Project số 3: Ảnh hưởng của tập thể dục với sức khỏe con người

1 .Import tất cả các library cần thiết

library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(janitor)
## 
## Attaching package: 'janitor'
## 
## The following objects are masked from 'package:stats':
## 
##     chisq.test, fisher.test
library(ggplot2)
library(lmPerm)
## Warning: package 'lmPerm' was built under R version 4.4.2
library(klaR)
## Warning: package 'klaR' was built under R version 4.4.2
## Loading required package: MASS
## 
## Attaching package: 'MASS'
## 
## The following object is masked from 'package:dplyr':
## 
##     select
library(MASS)
library(nnet)
library(caret) 
## Warning: package 'caret' was built under R version 4.4.2
## Loading required package: lattice
## 
## Attaching package: 'caret'
## 
## The following object is masked from 'package:purrr':
## 
##     lift
library(corrplot)
## corrplot 0.94 loaded
library(MASS)
library(rpart)
## Warning: package 'rpart' was built under R version 4.4.2
library(rpart.plot)
## Warning: package 'rpart.plot' was built under R version 4.4.2
library(xgboost)
## Warning: package 'xgboost' was built under R version 4.4.2
## 
## Attaching package: 'xgboost'
## 
## The following object is masked from 'package:dplyr':
## 
##     slice
library(Matrix)
## 
## Attaching package: 'Matrix'
## 
## The following objects are masked from 'package:tidyr':
## 
##     expand, pack, unpack
library(randomForest)
## Warning: package 'randomForest' was built under R version 4.4.2
## randomForest 4.7-1.2
## Type rfNews() to see new features/changes/bug fixes.
## 
## Attaching package: 'randomForest'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
## 
## The following object is masked from 'package:ggplot2':
## 
##     margin
library(ggcorrplot)
## Warning: package 'ggcorrplot' was built under R version 4.4.2
library(patchwork)
## Warning: package 'patchwork' was built under R version 4.4.2
## 
## Attaching package: 'patchwork'
## 
## The following object is masked from 'package:MASS':
## 
##     area
library(boot)
## 
## Attaching package: 'boot'
## 
## The following object is masked from 'package:lattice':
## 
##     melanoma
library(gridExtra)
## Warning: package 'gridExtra' was built under R version 4.4.2
## 
## Attaching package: 'gridExtra'
## 
## The following object is masked from 'package:randomForest':
## 
##     combine
## 
## The following object is masked from 'package:dplyr':
## 
##     combine
library(rpart.plot)
library(xgboost)
library(Matrix)
library(naivebayes)
## Warning: package 'naivebayes' was built under R version 4.4.2
## naivebayes 1.0.0 loaded
## For more information please visit: 
## https://majkamichal.github.io/naivebayes/
library(klaR)

2 .Khám phá và phân tích dữ liệu

2.1 .Tiền xử lí cơ bản

  • Đọc dữ liệu và làm sạch tên cột.
#Đọc dữ liệu
data = read.csv('https://raw.githubusercontent.com/vdnghia03/Dataset_AI_DS_DA_Repo/refs/heads/main/bodyPerformance.csv')
#Làm sạch tên cột
data = data |> clean_names()

glimpse(data)
## Rows: 13,393
## Columns: 12
## $ age                     <dbl> 27, 25, 31, 32, 28, 36, 42, 33, 54, 28, 42, 57…
## $ gender                  <chr> "M", "M", "M", "M", "M", "F", "F", "M", "M", "…
## $ height_cm               <dbl> 172.3, 165.0, 179.6, 174.5, 173.8, 165.4, 164.…
## $ weight_kg               <dbl> 75.24, 55.80, 78.00, 71.10, 67.70, 55.40, 63.7…
## $ body_fat                <dbl> 21.3, 15.7, 20.1, 18.4, 17.1, 22.0, 32.2, 36.9…
## $ diastolic               <dbl> 80, 77, 92, 76, 70, 64, 72, 84, 85, 81, 63, 69…
## $ systolic                <dbl> 130, 126, 152, 147, 127, 119, 135, 137, 165, 1…
## $ grip_force              <dbl> 54.9, 36.4, 44.8, 41.4, 43.5, 23.8, 22.7, 45.9…
## $ sit_and_bend_forward_cm <dbl> 18.4, 16.3, 12.0, 15.2, 27.1, 21.0, 0.8, 12.3,…
## $ sit_ups_counts          <dbl> 60, 53, 49, 53, 45, 27, 18, 42, 34, 55, 68, 0,…
## $ broad_jump_cm           <dbl> 217, 229, 181, 219, 217, 153, 146, 234, 148, 2…
## $ class                   <chr> "C", "A", "C", "B", "B", "B", "D", "B", "C", "…
  • Kiểm tra số lượng giá trị bị thiếu.
colSums(is.na(data))
##                     age                  gender               height_cm 
##                       0                       0                       0 
##               weight_kg                body_fat               diastolic 
##                       0                       0                       0 
##                systolic              grip_force sit_and_bend_forward_cm 
##                       0                       0                       0 
##          sit_ups_counts           broad_jump_cm                   class 
##                       0                       0                       0

Dữ liệu khá tốt, không bị vấn đề khuyết dữ liệu

Xử lí một số cột bị sai dữ liệu

Huyết áp tâm trương (diastolic) và tâm thu (systolic) không thể bằng 0 và tâm trương phải nhỏ hơn tâm thu

Xử lý dữ liệu cột diastolic và systolic = 0 và diastolic >= systolic

# Huyết áp tâm trương (diastolic) và tâm thu (systolic) không thể bằng 0 và
# tâm trương phải nhỏ hơn tâm thu
# Xử lý dữ liệu cột diastolic và systolic = 0 và diastolic >= systolic
data <- data %>%
  filter(diastolic != 0, systolic != 0, diastolic < systolic)

Xử lý dữ liệu cột grip_force và broad_jump_cm = 0

# Xử lý dữ liệu cột grip_force và broad_jump_cm = 0
data <- data %>%
  filter(grip_force != 0, broad_jump_cm != 0)
# Kiểm tra dữ liệu bị mất cân bằng
# Tính tổng số các class trong data
class_counts <- table(data$class)
print(class_counts)
## 
##    A    B    C    D 
## 3347 3345 3343 3340

Tính tổng số nam và nữ trong data

# Tính tổng số nam và nữ trong data
gender_counts <- table(data$gender)
print(gender_counts)
## 
##    F    M 
## 4916 8459

Qua đó ta thấy được dữ liệu bị mất cân bằng cột giới tính (gender) giữa nam và nữ (tỷ lệ 6.3 nam : 3.7 nữ) nhưng chưa tới mức quá lệch , nên ta vẫn sẽ dùng bộ dữ liệu ban đầu chứ không dùng các phương nhân bản dữ liệu

Chuyển các biến genderclass sang kiểu factor

data <- data %>%
  mutate(
    gender = as.factor(gender),
    class = as.factor(class)
  )

2.2 . Thêm cột, thêm các biến phân loại

Tạo biến mới là biến BMI rồi sau đó phân loại theo WHO

#Tạo biến mới tính BMI
data <- data %>%
  mutate(bmi = weight_kg / (height_cm / 100)^2)

# Phân loại BMI theo chuẩn WHO
data <- data %>%
  mutate(
    bmi_category = case_when(
      bmi < 18.5 ~ "Underweight",
      bmi >= 18.5 & bmi < 25 ~ "Normal",
      bmi >= 25 & bmi < 30 ~ "Overweight",
      bmi >= 30 ~ "Obese"
    )
  )

glimpse(data)
## Rows: 13,375
## Columns: 14
## $ age                     <dbl> 27, 25, 31, 32, 28, 36, 42, 33, 54, 28, 42, 57…
## $ gender                  <fct> M, M, M, M, M, F, F, M, M, M, M, F, F, M, M, F…
## $ height_cm               <dbl> 172.3, 165.0, 179.6, 174.5, 173.8, 165.4, 164.…
## $ weight_kg               <dbl> 75.24, 55.80, 78.00, 71.10, 67.70, 55.40, 63.7…
## $ body_fat                <dbl> 21.3, 15.7, 20.1, 18.4, 17.1, 22.0, 32.2, 36.9…
## $ diastolic               <dbl> 80, 77, 92, 76, 70, 64, 72, 84, 85, 81, 63, 69…
## $ systolic                <dbl> 130, 126, 152, 147, 127, 119, 135, 137, 165, 1…
## $ grip_force              <dbl> 54.9, 36.4, 44.8, 41.4, 43.5, 23.8, 22.7, 45.9…
## $ sit_and_bend_forward_cm <dbl> 18.4, 16.3, 12.0, 15.2, 27.1, 21.0, 0.8, 12.3,…
## $ sit_ups_counts          <dbl> 60, 53, 49, 53, 45, 27, 18, 42, 34, 55, 68, 0,…
## $ broad_jump_cm           <dbl> 217, 229, 181, 219, 217, 153, 146, 234, 148, 2…
## $ class                   <fct> C, A, C, B, B, B, D, B, C, B, A, D, C, C, C, A…
## $ bmi                     <dbl> 25.34418, 20.49587, 24.18143, 23.34956, 22.412…
## $ bmi_category            <chr> "Overweight", "Normal", "Normal", "Normal", "N…

Phân loại huyết áp dựa trên 2 biến diatolic và systolic

data <- data %>%
  mutate(
    blood_pressure_category = case_when(
      systolic < 90 & diastolic < 60 ~ "Hypotension", #Huyết áp thấp: systolic < 90 và diastolic < 60 mmHg
      systolic < 120 & diastolic < 80 ~ "Normal", # Huyết áp bình thường: systolic< 120 và diastolic < 80 mmHg
      systolic >= 120 & systolic < 130 & diastolic < 80 ~ "Elevated", #Huyết áp cao
      systolic >= 130 & systolic < 140 | diastolic >= 80 & diastolic < 90 ~ "Hypertension Stage 1", #Huyết áp cao giai đoạn 1
      systolic >= 140 | diastolic >= 90 ~ "Hypertension Stage 2", #Huyết ́ pao giai đoạn 2
      systolic > 180 | diastolic > 120 ~ "Hypertensive Crisis" #Huyết áp ở mức nguy hiểm
    )
  )

Chuyển các biến bmi_categoryblood_pressure_category sang kiểu factor

data <- data %>%
  mutate(
    bmi_category = as.factor(bmi_category),
    blood_pressure_category = as.factor(blood_pressure_category)
  )

Tạo biến mới là biến finess score và biến hiệu số huyết áp

  • Biến finess score: Dựa trên các yếu tố thể lực như sit_ups_counts, grip_force và broad_jump_cm với các hệ số để tạo nên

  • Biến pulse_pressure: Vì các biến như systolic và diastolic chỉ thể hiện tại một thời điểm nhất định, nên ta tạo ra biến hiệu số huyết áp bp_difference để nhìn rõ hơn.

data <- data %>%
  mutate(
    # Fitness Score dựa trên các yếu tố thể lực
    fitness_score = sit_ups_counts * 0.4 + grip_force * 0.3 + broad_jump_cm * 0.3,
    # Hiệu số huyết áp
    pulse_pressure = systolic - diastolic,
  )
################
glimpse(data)
## Rows: 13,375
## Columns: 17
## $ age                     <dbl> 27, 25, 31, 32, 28, 36, 42, 33, 54, 28, 42, 57…
## $ gender                  <fct> M, M, M, M, M, F, F, M, M, M, M, F, F, M, M, F…
## $ height_cm               <dbl> 172.3, 165.0, 179.6, 174.5, 173.8, 165.4, 164.…
## $ weight_kg               <dbl> 75.24, 55.80, 78.00, 71.10, 67.70, 55.40, 63.7…
## $ body_fat                <dbl> 21.3, 15.7, 20.1, 18.4, 17.1, 22.0, 32.2, 36.9…
## $ diastolic               <dbl> 80, 77, 92, 76, 70, 64, 72, 84, 85, 81, 63, 69…
## $ systolic                <dbl> 130, 126, 152, 147, 127, 119, 135, 137, 165, 1…
## $ grip_force              <dbl> 54.9, 36.4, 44.8, 41.4, 43.5, 23.8, 22.7, 45.9…
## $ sit_and_bend_forward_cm <dbl> 18.4, 16.3, 12.0, 15.2, 27.1, 21.0, 0.8, 12.3,…
## $ sit_ups_counts          <dbl> 60, 53, 49, 53, 45, 27, 18, 42, 34, 55, 68, 0,…
## $ broad_jump_cm           <dbl> 217, 229, 181, 219, 217, 153, 146, 234, 148, 2…
## $ class                   <fct> C, A, C, B, B, B, D, B, C, B, A, D, C, C, C, A…
## $ bmi                     <dbl> 25.34418, 20.49587, 24.18143, 23.34956, 22.412…
## $ bmi_category            <fct> Overweight, Normal, Normal, Normal, Normal, No…
## $ blood_pressure_category <fct> Hypertension Stage 1, Elevated, Hypertension S…
## $ fitness_score           <dbl> 105.57, 100.82, 87.34, 99.32, 96.15, 63.84, 57…
## $ pulse_pressure          <dbl> 50, 49, 60, 71, 57, 55, 63, 53, 80, 75, 47, 37…

Mô tả và biểu diễn tổng hợp dữ liệu

# Tóm tắt dữ liệu cơ bản
summary(data)
##       age        gender     height_cm       weight_kg         body_fat    
##  Min.   :21.00   F:4916   Min.   :125.0   Min.   : 26.30   Min.   : 3.00  
##  1st Qu.:25.00   M:8459   1st Qu.:162.4   1st Qu.: 58.20   1st Qu.:18.00  
##  Median :32.00            Median :169.2   Median : 67.44   Median :22.80  
##  Mean   :36.77            Mean   :168.6   Mean   : 67.45   Mean   :23.24  
##  3rd Qu.:48.00            3rd Qu.:174.8   3rd Qu.: 75.30   3rd Qu.:28.00  
##  Max.   :64.00            Max.   :193.8   Max.   :138.10   Max.   :78.40  
##    diastolic        systolic       grip_force    sit_and_bend_forward_cm
##  Min.   :  6.0   Min.   : 77.0   Min.   : 1.60   Min.   :-25.00         
##  1st Qu.: 71.0   1st Qu.:120.0   1st Qu.:27.50   1st Qu.: 10.90         
##  Median : 79.0   Median :130.0   Median :37.90   Median : 16.20         
##  Mean   : 78.8   Mean   :130.3   Mean   :36.98   Mean   : 15.21         
##  3rd Qu.: 86.0   3rd Qu.:141.0   3rd Qu.:45.20   3rd Qu.: 20.75         
##  Max.   :126.0   Max.   :201.0   Max.   :70.50   Max.   :213.00         
##  sit_ups_counts  broad_jump_cm   class         bmi             bmi_category 
##  Min.   : 0.00   Min.   : 20.0   A:3347   Min.   :11.10   Normal     :9233  
##  1st Qu.:30.00   1st Qu.:162.0   B:3345   1st Qu.:21.61   Obese      : 342  
##  Median :41.00   Median :193.0   C:3343   Median :23.46   Overweight :3463  
##  Mean   :39.78   Mean   :190.3   D:3340   Mean   :23.61   Underweight: 337  
##  3rd Qu.:50.00   3rd Qu.:221.0            3rd Qu.:25.34                     
##  Max.   :80.00   Max.   :303.0            Max.   :42.91                     
##          blood_pressure_category fitness_score    pulse_pressure  
##  Elevated            :2089       Min.   : 11.38   Min.   :  6.00  
##  Hypertension Stage 1:5929       1st Qu.: 71.08   1st Qu.: 44.00  
##  Hypertension Stage 2:2465       Median : 85.79   Median : 51.00  
##  Hypotension         :   6       Mean   : 84.09   Mean   : 51.46  
##  Normal              :2886       3rd Qu.: 99.39   3rd Qu.: 59.00  
##                                  Max.   :133.30   Max.   :139.00

3 .PHẦN 1 : BẢNG TỔNG HỢP VÀ MỘT VÀI BIỂU ĐỒ TRỰC QUAN DỮ LIỆU ĐỒNG THỜI ÁP DỤNG A/B TESTINGS ĐỂ KIỂM ĐỊNH CÁC GIẢ THUYẾT CẦN THIẾT

3.1 .Tóm tắt dữ liệu theo nhóm phân loại gender và class

a. Theo gender

# Tính toán thống kê cho dữ liệu theo giới tính
pivot_gender <- data %>%
  group_by(gender) %>%
  summarise(
    mean_height = mean(height_cm, na.rm = TRUE),
    sd_height = sd(height_cm, na.rm = TRUE),
    mean_weight = mean(weight_kg, na.rm = TRUE),
    sd_weight = sd(weight_kg, na.rm = TRUE),
    mean_body_fat = mean(body_fat, na.rm = TRUE),
    sd_body_fat = sd(body_fat, na.rm = TRUE),
    mean_diastolic = mean(diastolic, na.rm = TRUE),
    sd_diastolic = sd(diastolic, na.rm = TRUE),
    mean_systolic = mean(systolic, na.rm = TRUE),
    sd_systolic = sd(systolic, na.rm = TRUE),
    mean_grip_force = mean(grip_force, na.rm = TRUE),
    sd_grip_force = sd(grip_force, na.rm = TRUE),
    mean_sit_and_bend_forward = mean(sit_and_bend_forward_cm, na.rm = TRUE),
    sd_sit_and_bend_forward = sd(sit_and_bend_forward_cm, na.rm = TRUE),
    mean_sit_ups = mean(sit_ups_counts, na.rm = TRUE),
    sd_sit_ups = sd(sit_ups_counts, na.rm = TRUE),
    mean_broad_jump = mean(broad_jump_cm, na.rm = TRUE),
    sd_broad_jump = sd(broad_jump_cm, na.rm = TRUE),
    mean_bmi = mean(bmi, na.rm = TRUE),
    sd_bmi = sd(bmi, na.rm = TRUE),
    mean_fitness_score = mean(fitness_score, na.rm = TRUE),
    sd_fitness_score = sd(fitness_score, na.rm = TRUE)
  )

glimpse(data)
## Rows: 13,375
## Columns: 17
## $ age                     <dbl> 27, 25, 31, 32, 28, 36, 42, 33, 54, 28, 42, 57…
## $ gender                  <fct> M, M, M, M, M, F, F, M, M, M, M, F, F, M, M, F…
## $ height_cm               <dbl> 172.3, 165.0, 179.6, 174.5, 173.8, 165.4, 164.…
## $ weight_kg               <dbl> 75.24, 55.80, 78.00, 71.10, 67.70, 55.40, 63.7…
## $ body_fat                <dbl> 21.3, 15.7, 20.1, 18.4, 17.1, 22.0, 32.2, 36.9…
## $ diastolic               <dbl> 80, 77, 92, 76, 70, 64, 72, 84, 85, 81, 63, 69…
## $ systolic                <dbl> 130, 126, 152, 147, 127, 119, 135, 137, 165, 1…
## $ grip_force              <dbl> 54.9, 36.4, 44.8, 41.4, 43.5, 23.8, 22.7, 45.9…
## $ sit_and_bend_forward_cm <dbl> 18.4, 16.3, 12.0, 15.2, 27.1, 21.0, 0.8, 12.3,…
## $ sit_ups_counts          <dbl> 60, 53, 49, 53, 45, 27, 18, 42, 34, 55, 68, 0,…
## $ broad_jump_cm           <dbl> 217, 229, 181, 219, 217, 153, 146, 234, 148, 2…
## $ class                   <fct> C, A, C, B, B, B, D, B, C, B, A, D, C, C, C, A…
## $ bmi                     <dbl> 25.34418, 20.49587, 24.18143, 23.34956, 22.412…
## $ bmi_category            <fct> Overweight, Normal, Normal, Normal, Normal, No…
## $ blood_pressure_category <fct> Hypertension Stage 1, Elevated, Hypertension S…
## $ fitness_score           <dbl> 105.57, 100.82, 87.34, 99.32, 96.15, 63.84, 57…
## $ pulse_pressure          <dbl> 50, 49, 60, 71, 57, 55, 63, 53, 80, 75, 47, 37…

b. Phân loại class

# Tóm tắt dữ liệu theo nhóm phân loại class

pivot_class <- data %>%
  group_by(class) %>%
  summarise(
    mean_height = mean(height_cm, na.rm = TRUE),
    sd_height = sd(height_cm, na.rm = TRUE),
    mean_weight = mean(weight_kg, na.rm = TRUE),
    sd_weight = sd(weight_kg, na.rm = TRUE),
    mean_body_fat = mean(body_fat, na.rm = TRUE),
    sd_body_fat = sd(body_fat, na.rm = TRUE),
    mean_diastolic = mean(diastolic, na.rm = TRUE),
    sd_diastolic = sd(diastolic, na.rm = TRUE),
    mean_systolic = mean(systolic, na.rm = TRUE),
    sd_systolic = sd(systolic, na.rm = TRUE),
    mean_grip_force = mean(grip_force, na.rm = TRUE),
    sd_grip_force = sd(grip_force, na.rm = TRUE),
    mean_sit_and_bend_forward = mean(sit_and_bend_forward_cm, na.rm = TRUE),
    sd_sit_and_bend_forward = sd(sit_and_bend_forward_cm, na.rm = TRUE),
    mean_sit_ups = mean(sit_ups_counts, na.rm = TRUE),
    sd_sit_ups = sd(sit_ups_counts, na.rm = TRUE),
    mean_broad_jump = mean(broad_jump_cm, na.rm = TRUE),
    sd_broad_jump = sd(broad_jump_cm, na.rm = TRUE),
    mean_bmi = mean(bmi, na.rm = TRUE),
    sd_bmi = sd(bmi, na.rm = TRUE),
    mean_fitness_score = mean(fitness_score, na.rm = TRUE),
    sd_fitness_score = sd(fitness_score, na.rm = TRUE)
  )

glimpse(data)
## Rows: 13,375
## Columns: 17
## $ age                     <dbl> 27, 25, 31, 32, 28, 36, 42, 33, 54, 28, 42, 57…
## $ gender                  <fct> M, M, M, M, M, F, F, M, M, M, M, F, F, M, M, F…
## $ height_cm               <dbl> 172.3, 165.0, 179.6, 174.5, 173.8, 165.4, 164.…
## $ weight_kg               <dbl> 75.24, 55.80, 78.00, 71.10, 67.70, 55.40, 63.7…
## $ body_fat                <dbl> 21.3, 15.7, 20.1, 18.4, 17.1, 22.0, 32.2, 36.9…
## $ diastolic               <dbl> 80, 77, 92, 76, 70, 64, 72, 84, 85, 81, 63, 69…
## $ systolic                <dbl> 130, 126, 152, 147, 127, 119, 135, 137, 165, 1…
## $ grip_force              <dbl> 54.9, 36.4, 44.8, 41.4, 43.5, 23.8, 22.7, 45.9…
## $ sit_and_bend_forward_cm <dbl> 18.4, 16.3, 12.0, 15.2, 27.1, 21.0, 0.8, 12.3,…
## $ sit_ups_counts          <dbl> 60, 53, 49, 53, 45, 27, 18, 42, 34, 55, 68, 0,…
## $ broad_jump_cm           <dbl> 217, 229, 181, 219, 217, 153, 146, 234, 148, 2…
## $ class                   <fct> C, A, C, B, B, B, D, B, C, B, A, D, C, C, C, A…
## $ bmi                     <dbl> 25.34418, 20.49587, 24.18143, 23.34956, 22.412…
## $ bmi_category            <fct> Overweight, Normal, Normal, Normal, Normal, No…
## $ blood_pressure_category <fct> Hypertension Stage 1, Elevated, Hypertension S…
## $ fitness_score           <dbl> 105.57, 100.82, 87.34, 99.32, 96.15, 63.84, 57…
## $ pulse_pressure          <dbl> 50, 49, 60, 71, 57, 55, 63, 53, 80, 75, 47, 37…

3.2 . Tạo bảng tổng hợp các biến và xem thống kê

Tạo bảng tổng hợp các biến

#Tạo bảng tổng hợp các biến
df_sum <- data %>%
  summarise(across(c(age, height_cm, weight_kg, body_fat, diastolic, systolic, grip_force, 
                     sit_and_bend_forward_cm, sit_ups_counts, broad_jump_cm ,bmi, fitness_score ,pulse_pressure),
                   list(gtnn = ~min(., na.rm = TRUE),
                        gtln = ~max(., na.rm = TRUE),
                        tv = ~median(., na.rm = TRUE),
                        tb = ~mean(., na.rm = TRUE),
                        dlc = ~sd(., na.rm = TRUE))))
df_sum
##   age_gtnn age_gtln age_tv   age_tb age_dlc height_cm_gtnn height_cm_gtln
## 1       21       64     32 36.76972  13.623            125          193.8
##   height_cm_tv height_cm_tb height_cm_dlc weight_kg_gtnn weight_kg_gtln
## 1        169.2     168.5643       8.42634           26.3          138.1
##   weight_kg_tv weight_kg_tb weight_kg_dlc body_fat_gtnn body_fat_gtln
## 1        67.44     67.44829       11.9471             3          78.4
##   body_fat_tv body_fat_tb body_fat_dlc diastolic_gtnn diastolic_gtln
## 1        22.8    23.23516     7.251582              6            126
##   diastolic_tv diastolic_tb diastolic_dlc systolic_gtnn systolic_gtln
## 1           79     78.80013       10.6945            77           201
##   systolic_tv systolic_tb systolic_dlc grip_force_gtnn grip_force_gtln
## 1         130    130.2641     14.61555             1.6            70.5
##   grip_force_tv grip_force_tb grip_force_dlc sit_and_bend_forward_cm_gtnn
## 1          37.9       36.9791        10.6052                          -25
##   sit_and_bend_forward_cm_gtln sit_and_bend_forward_cm_tv
## 1                          213                       16.2
##   sit_and_bend_forward_cm_tb sit_and_bend_forward_cm_dlc sit_ups_counts_gtnn
## 1                   15.20959                    8.457525                   0
##   sit_ups_counts_gtln sit_ups_counts_tv sit_ups_counts_tb sit_ups_counts_dlc
## 1                  80                41          39.78437           14.25838
##   broad_jump_cm_gtnn broad_jump_cm_gtln broad_jump_cm_tv broad_jump_cm_tb
## 1                 20                303              193         190.2683
##   broad_jump_cm_dlc bmi_gtnn bmi_gtln   bmi_tv   bmi_tb  bmi_dlc
## 1          39.54228 11.10398 42.90651 23.46098 23.60516 2.939601
##   fitness_score_gtnn fitness_score_gtln fitness_score_tv fitness_score_tb
## 1              11.38              133.3            85.79         84.08796
##   fitness_score_dlc pulse_pressure_gtnn pulse_pressure_gtln pulse_pressure_tv
## 1          19.04238                   6                 139                51
##   pulse_pressure_tb pulse_pressure_dlc
## 1          51.46395           10.76123

Chuyển đổi và tách tên biến, giá trị thống kê

df_stats_tidy <- df_sum %>%
  pivot_longer(cols = everything(), names_to = "ten", values_to = "gt") %>% 
  separate(ten, into = c("bien", "stat"), sep = "_(?=[^_]+$)") %>% 
  pivot_wider(names_from = stat, values_from = gt)
# Kết quả
print(df_stats_tidy)
## # A tibble: 13 × 6
##    bien                     gtnn  gtln    tv    tb   dlc
##    <chr>                   <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 age                      21    64    32    36.8 13.6 
##  2 height_cm               125   194.  169.  169.   8.43
##  3 weight_kg                26.3 138.   67.4  67.4 11.9 
##  4 body_fat                  3    78.4  22.8  23.2  7.25
##  5 diastolic                 6   126    79    78.8 10.7 
##  6 systolic                 77   201   130   130.  14.6 
##  7 grip_force                1.6  70.5  37.9  37.0 10.6 
##  8 sit_and_bend_forward_cm -25   213    16.2  15.2  8.46
##  9 sit_ups_counts            0    80    41    39.8 14.3 
## 10 broad_jump_cm            20   303   193   190.  39.5 
## 11 bmi                      11.1  42.9  23.5  23.6  2.94
## 12 fitness_score            11.4 133.   85.8  84.1 19.0 
## 13 pulse_pressure            6   139    51    51.5 10.8
print(pivot_gender)
## # A tibble: 2 × 23
##   gender mean_height sd_height mean_weight sd_weight mean_body_fat sd_body_fat
##   <fct>        <dbl>     <dbl>       <dbl>     <dbl>         <dbl>       <dbl>
## 1 F             160.      5.65        56.9      7.63          28.5        6.22
## 2 M             173.      5.81        73.6      9.47          20.2        5.95
## # ℹ 16 more variables: mean_diastolic <dbl>, sd_diastolic <dbl>,
## #   mean_systolic <dbl>, sd_systolic <dbl>, mean_grip_force <dbl>,
## #   sd_grip_force <dbl>, mean_sit_and_bend_forward <dbl>,
## #   sd_sit_and_bend_forward <dbl>, mean_sit_ups <dbl>, sd_sit_ups <dbl>,
## #   mean_broad_jump <dbl>, sd_broad_jump <dbl>, mean_bmi <dbl>, sd_bmi <dbl>,
## #   mean_fitness_score <dbl>, sd_fitness_score <dbl>
print(pivot_class)
## # A tibble: 4 × 23
##   class mean_height sd_height mean_weight sd_weight mean_body_fat sd_body_fat
##   <fct>       <dbl>     <dbl>       <dbl>     <dbl>         <dbl>       <dbl>
## 1 A            168.      7.84        64.4      10.6          20.5        6.44
## 2 B            169.      8.13        66.6      10.9          22.0        6.65
## 3 C            169.      8.52        66.8      10.9          22.6        6.27
## 4 D            169.      9.11        72.0      13.9          27.7        7.50
## # ℹ 16 more variables: mean_diastolic <dbl>, sd_diastolic <dbl>,
## #   mean_systolic <dbl>, sd_systolic <dbl>, mean_grip_force <dbl>,
## #   sd_grip_force <dbl>, mean_sit_and_bend_forward <dbl>,
## #   sd_sit_and_bend_forward <dbl>, mean_sit_ups <dbl>, sd_sit_ups <dbl>,
## #   mean_broad_jump <dbl>, sd_broad_jump <dbl>, mean_bmi <dbl>, sd_bmi <dbl>,
## #   mean_fitness_score <dbl>, sd_fitness_score <dbl>

Phát hiện giá trị sit_and_bend_forward_cm = 213 là không hợp lí , vì giá trị cao nhất của chỉ số này chỉ thường ở mức >46cm và <50cm

Loại bỏ các giá trị > 50cm trong cột sit_and_bend_forward_cm

data <- data[data$sit_and_bend_forward_cm <= 50, ]
# xem lại kết quả : 
df_sum <- data %>%
  summarise(across(c(age, height_cm, weight_kg, body_fat, diastolic, systolic, grip_force, 
                     sit_and_bend_forward_cm, sit_ups_counts, broad_jump_cm ,bmi, fitness_score ,pulse_pressure),
                   list(gtnn = ~min(., na.rm = TRUE),
                        gtln = ~max(., na.rm = TRUE),
                        tv = ~median(., na.rm = TRUE),
                        tb = ~mean(., na.rm = TRUE),
                        dlc = ~sd(., na.rm = TRUE))))

Chuyển đổi và tách tên biến, giá trị thống kê

df_stats_tidy <- df_sum %>%
  pivot_longer(cols = everything(), names_to = "ten", values_to = "gt") %>% 
  separate(ten, into = c("bien", "stat"), sep = "_(?=[^_]+$)") %>% 
  pivot_wider(names_from = stat, values_from = gt)

print(df_stats_tidy)
## # A tibble: 13 × 6
##    bien                     gtnn  gtln    tv    tb   dlc
##    <chr>                   <dbl> <dbl> <dbl> <dbl> <dbl>
##  1 age                      21    64    32    36.8 13.6 
##  2 height_cm               125   194.  169.  169.   8.43
##  3 weight_kg                26.3 138.   67.4  67.4 11.9 
##  4 body_fat                  3    78.4  22.8  23.2  7.25
##  5 diastolic                 6   126    79    78.8 10.7 
##  6 systolic                 77   201   130   130.  14.6 
##  7 grip_force                1.6  70.5  37.9  37.0 10.6 
##  8 sit_and_bend_forward_cm -25    42    16.2  15.2  8.15
##  9 sit_ups_counts            0    80    41    39.8 14.3 
## 10 broad_jump_cm            20   303   193   190.  39.5 
## 11 bmi                      11.1  42.9  23.5  23.6  2.94
## 12 fitness_score            11.4 133.   85.8  84.1 19.0 
## 13 pulse_pressure            6   139    51    51.5 10.8

3.3 .Vẽ một vài biểu đồ trực quan nhằm khai thác một số khía cạnh của dữ liệu

Tạo histogram cho từng biến

numeric_cols <- dplyr::select(data, age, height_cm, weight_kg, body_fat, 
                       diastolic, systolic, grip_force, sit_and_bend_forward_cm,
                       sit_ups_counts, broad_jump_cm, bmi, fitness_score, pulse_pressure)

# Kiểm tra phân phối cho từng biến bằng histogram
plots <- list()
for (col in names(numeric_cols)) {
  plot <- ggplot(data, aes_string(x = col)) +
    geom_histogram(binwidth = 2, fill = "blue", color = "white") +
    labs(
      x = col,
      y = NULL  # Xóa tiêu đề trục Oy
    ) +
    theme_minimal() +
    theme(
      axis.title.y = element_blank(),  # Xóa nhãn trục Oy
      axis.ticks.y = element_blank()   # Xóa dấu tick trên trục Oy
    )
  plots[[col]] <- plot
}
## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
# Ghép tất cả các biểu đồ histogram thành một hình
combined_plot <- wrap_plots(plots, ncol = 4)  # Số cột, ví dụ 3
print(combined_plot)

3.4 .So sánh giữa các nhóm

3.4.1 .Theo giới tính

1. grip_force dựa trên class và gender

ggplot(data, aes(x = class, y = grip_force, fill = gender)) +
  geom_bar(stat = "summary", fun = "mean", position = "dodge") +
  labs(title = "Lực kẹp trung bình theo nhóm hiệu suất và giới tính", 
       x = "Nhóm hiệu suất", y = "Lực kẹp (kg)") +
  theme_minimal()

Nhận xét : Chỉ số lực kẹp ở nam giới nhìn chung tốt hơn nữ giới . Chỉ số này tăng đều theo phân lớp hiệu quả (class)

Ta đi kiểm định nhận xét trên

# Hàm permutation test với 2 nhóm
perm_fun <- function(x, nA, nB, R) {
  n <- nA + nB
  mean_diff <- numeric(R)
  for (i in 1:R){
    idx_a <- sample(x = 1:n, size = nA)
    idx_b <- setdiff(x = 1:n, y = idx_a)
    mean_diff[i] <- mean(x[idx_a]) - mean(x[idx_b])
  }
  return(mean_diff)
}
class_counts <- table(data$class)
print(class_counts)
## 
##    A    B    C    D 
## 3346 3344 3343 3340
# Tính tổng số nam và nữ trong data
gender_counts <- table(data$gender)
print(gender_counts)
## 
##    F    M 
## 4916 8457

Để kiểm tra xem nhận xét mang ý nghĩa thống kê hay chỉ là kết quả của sự ngẫu nhiên , ta thực hiện các kiểm định sau : H0: Không có sự khác biệt đáng kể về lực kẹp trung bình giữa nam và nữ H1: Lực kẹp trung bình ở nữ yếu hơn nam

diff_h1 <- perm_fun(x = data$grip_force, nA = 4916, nB = 8457, R = 1000)
h1_mean_a <- mean(data$grip_force[data$gender == "F"])
h1_mean_b <- mean(data$grip_force[data$gender == "M"])
mean(diff_h1 < (h1_mean_a - h1_mean_b))
## [1] 0

P-value = 0 Với mức ý nghĩa alpha = 0.05 , p-value < alpha nên chấp nhận H1 : Lực kẹp trung bình ở nữ yếu hơn nam

Kế tiếp ta có

H0: Không có sự khác biệt đáng kể về lực kẹp trung bình giữa các nhóm hiệu suất. H1: Có ít nhất một sự khác biệt đáng kể về lực kẹp trung bình giữa các nhóm hiệu suất.

h1 <- aovp(formula = grip_force ~ class, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h1)
## Component 1 :
##                Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## class           3    28475    9491.7 5000 < 2.2e-16 ***
## Residuals   13369  1475698     110.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value = 2.2e-16 Với mức ý nghĩa alpha = 0.05 , p-val < alpha nên ta chấp nhận H1 : Có ít nhất một sự khác biệt đáng kể về lực kẹp trung bình giữa các nhóm hiệu suất.
Vậy kết luận nhận xét ban đầu là chính xác và mang ý nghĩa thống kê

2. heigh_cm dựa trên bmi_category và gender

ggplot(data, aes(x = bmi_category, y = height_cm, fill = gender)) +
  geom_bar(stat = "summary", fun = "mean", position = "dodge") +
  labs(title = "Chiều cao trung bình theo phân loại BMI và giới tính", 
       x = "Phân loại BMI", y = "Chiều cao (cm)") +
  theme_minimal()

Nhận xét : Chiều cao trung bình ở nam giới trong cả 4 phân nhóm BMI là gần như ngang nhau và tốt hơn nữ giới , trong đó nhóm Underweight có chỉ số chiều cao kém hơn một chút . Ở nữ giới thì 2 nhóm Obese và Overweight thì lại có chiều cao trung bình thấp hơn so với 2 nhóm Normal và Underweight
Để kiểm tra xem nhận xét mang ý nghĩa thống kê hay chỉ là kết quả của sự ngẫu nhiên , ta thực hiện các kiểm định sau :

H0: Không có sự khác biệt đáng kể về chiều cao trung bình giữa nam và nữ
H1: chiều cao trung bình ở nữ thấp hơn nam

diff_h2 <- perm_fun(x = data$height_cm, nA = 4916, nB = 8457, R = 1000)
h2_mean_a <- mean(data$height_cm[data$gender == "F"])
h2_mean_b <- mean(data$height_cm[data$gender == "M"])
mean(diff_h2 < (h2_mean_a - h2_mean_b))
## [1] 0

P-value = 0 Với mức ý nghĩa alpha = 0.05 , p-value < alpha nên chấp nhận H1 : chiều cao trung bình ở nữ thấp hơn nam
Kế tiếp ta có :
H0: Không có sự khác biệt về chiều cao trung bình giữa các nhóm BMI ở nam giới.
H1: Có ít nhất một nhóm BMI có chiều cao trung bình khác biệt ở nam giới.

# Lọc dữ liệu chỉ lấy nam giới
data_male <- subset(data, gender == 'M')
h2 <- aovp(formula = height_cm ~ bmi_category, data = data_male, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h2)
## Component 1 :
##                Df R Sum Sq R Mean Sq Iter Pr(Prob)
## bmi_category    3      247    82.317 2464   0.6814
## Residuals    8453   285277    33.749

P-val = 0.8924
Với mức ý nghĩa alpha = 0.05 ,p-val > alpha nên chấp nhận H0 : Không có sự khác biệt về chiều cao trung bình giữa các nhóm BMI ở nam giới.

Kế tiếp ta có : H0: Không có sự khác biệt về chiều cao trung bình giữa nhóm {Obese, Overweight} và nhóm {Normal, Underweight}.
H1: Chiều cao trung bình giữa nhóm {Obese, Overweight} thấp hơn nhóm {Normal, Underweight}. Gộp nhóm {Obese, Overweight} và {Normal, Underweight}

data_female <- subset(data, gender == 'F')
# Gộp nhóm {Obese, Overweight} và {Normal, Underweight}

data_female$group <- ifelse(data_female$bmi_category %in% c("Obese", "Overweight"), 
                            "High BMI", "Low BMI")

print(table(data_female$group))
## 
## High BMI  Low BMI 
##      622     4294
summary(data_female$height_cm)
##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##   125.0   156.6   160.5   160.5   164.3   179.0
h2f <- aovp(formula = height_cm ~ group, data = data_female, perm = "Prob")
## [1] "Settings:  unique SS "
print(summary(h2f))
## Component 1 :
##               Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## group1         1     2732   2732.42 5000 < 2.2e-16 ***
## Residuals   4914   154216     31.38                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-val = 2.2e-16
Với mức ý nghĩa alpha =0.05 , p-val < alpha nên ta chấp nhận H1: Chiều cao trung bình giữa nhóm {Obese, Overweight} thấp hơn nhóm {Normal, Underweight}.
Vậy kết luận nhận xét ban đầu là chính xác và mang ý nghĩa thống kê

3.body_fat dựa trên class và gender

ggplot(data, aes(x = class, y = body_fat, fill = gender)) +
  geom_bar(stat = "summary", fun = "mean", position = "dodge") +
  labs(title = "Tỷ lệ mỡ cơ thể trung bình theo nhóm hiệu suất và giới tính", 
       x = "Nhóm hiệu suất", y = "Tỷ lệ mỡ (%)") +
  theme_minimal()

Nhận xét : % Mỡ (body-fat) thấp nhất ở class A , sau đó tăng dần đều ( ở cả 2 giới ) tới class D là cao nhất . % mỡ ở nam giới có xu hướng thấp hơn nữ giới

Để kiểm tra xem nhận xét mang ý nghĩa thống kê hay chỉ là kết quả của sự ngẫu nhiên , ta thực hiện các kiểm định:
H0: Không có sự khác biệt về chiều tỉ lệ mỡ trung bình ở nam giới và nữ giới.
H1: Tỉ lệ mỡ trung bình ở nam giới thấp hơn nữ giới.

diff_h3 <- perm_fun(x = data$body_fat, nA = 4916, nB = 8457, R = 1000)
h3_mean_a <- mean(data$body_fat[data$gender == "M"])
h3_mean_b <- mean(data$body_fat[data$gender == "F"])
mean(diff_h3 < (h3_mean_a - h3_mean_b))
## [1] 0

P-value = 0 Với mức ý nghĩa alpha = 0.05 , p-value < alpha nên chấp nhận H1 : Tỉ lệ mỡ cơ thể trung bình ở nam giới thấp hơn nữ giới Kế tiếp ta có :
H0: Không có sự khác biệt về tỉ lệ mỡ cơ thể trung bình ở các nhóm hiệu suất.
H1: Có ít nhất một nhóm hiệu suất có thỉ lệ mỡ cơ thể trung bình khác các nhóm còn lại.

h3 <- aovp(formula = body_fat ~ class, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h3)
## Component 1 :
##                Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## class           3    97546     32515 5000 < 2.2e-16 ***
## Residuals   13369   605701        45                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-val = 2.2e-16
Với mức ý nghĩa alpha = 0.05 , p-value < alpha nên chấp nhận H1: Có ít nhất một nhóm hiệu suất có tỉ lệ mỡ cơ thể trung bình khác các nhóm còn lại.
Vậy kết luận nhận xét ban đầu là chính xác và mang ý nghĩa thống kê

4. Vẽ biểu đồ trung bình fitness_score theo class và giới tính

ggplot(data, aes(x = class, y = fitness_score, fill = gender)) +
  geom_bar(stat = "summary", fun = "mean", position = "dodge") +
  labs(title = "Điểm sức khỏe trung bình theo nhóm hiệu suất và giới tính", 
       x = "Nhóm hiệu suất", y = "Điểm sức khỏe") +
  theme_minimal()

Nhận xét : Chỉ số điểm fitness score tốt nhất ở class A , sau đó giảm dần (ở cả 2 giới).Fitness score ở nam giới cao hơn nữ giới
Để kiểm tra xem nhận xét mang ý nghĩa thống kê hay chỉ là kết quả của sự ngẫu nhiên , ta thực hiện các kiểm định:
H0: Không có sự khác biệt về fitness score trung bình ở nam giới và nữ giới.
H1: fitness score trung bình ở nữ giới thấp hơn nam giới.

diff_h4 <- perm_fun(x = data$fitness_score, nA = 4916, nB = 8457, R = 1000)
h4_mean_a <- mean(data$fitness_score[data$gender == "M"])
h4_mean_b <- mean(data$fitness_score[data$gender == "F"])
mean(diff_h4 < (h4_mean_b - h4_mean_a))
## [1] 0

P-value = 0
Với mức ý nghĩa alpha = 0.05 , p-value < alpha nên chấp nhận H1 : fitness score trung bình ở nữ giới thấp hơn nam giới
Kế tiếp ta có :
H0: Không có sự khác biệt về fitness score trung bình ở các nhóm hiệu suất.
H1: Có ít nhất một nhóm hiệu suất có fitness score trung bình khác các nhóm còn lại.

h4 <- aovp(formula = fitness_score ~ class, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h4)
## Component 1 :
##                Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## class           3   515500    171833 5000 < 2.2e-16 ***
## Residuals   13369  4333453       324                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-val = 2.2e-16 Với mức ý nghĩa alpha = 0.05 , p-value < alpha nên chấp nhận H1: Có ít nhất một nhóm hiệu suất có fitness score trung bình khác các nhóm còn lại. Vậy kết luận nhận xét ban đầu là chính xác và mang ý nghĩa thống kê

5. Vẽ biểu đồ trung bình pulse_pressure theo class và giới tính

ggplot(data, aes(x = class, y = pulse_pressure, fill = gender)) +
  geom_bar(stat = "summary", fun = "mean", position = "dodge") +
  labs(title = "Sự khác biệt chỉ số Pulse_pressure(áp lực mạch đập) trung bình theo nhóm hiệu suất và giới tính", 
       x = "Nhóm hiệu suất", y = "Pulse pressure (mmHg)") +
  theme_minimal()

Nhận xét : Chỉ số pulse pressure trung bình ở nam giới giảm dần từ class A về class D và cao hơn nữ giới . Trong khi ở nữ giới lại không có xu hướng tăng giảm rõ ràng với class C thấp nhất , B cao nhất , A và D gần ngang nhau Để kiểm tra nhận xét trên có mang ý nghĩa thống kê không , ta thực hiện các kiểm định: H0 : Không có sự khác biệt về pulse pressure trung bình ở 2 giới tính H1 : Pulse pressure trung bình ở nữ giới thấp hơn nam giới

diff_h5 <- perm_fun(x = data$fitness_score, nA = 4916, nB = 8457, R = 1000)
h5_mean_a <- mean(data$pulse_pressure[data$gender == "F"])
h5_mean_b <- mean(data$pulse_pressure[data$gender == "M"])
mean(diff_h5 < (h5_mean_a - h5_mean_b))
## [1] 0
#p-val = 0 
#Với mức ý nghĩa alpha =0.05 , p-val < alpha nên chấp nhận H1 
#Kế tiếp ta có : 
#H0 : Không có sự khác biệt về pulse pressure trung bình giữa các nhóm hiệu suất của nam giới 
#H1 : Có ít nhất 1 sự khác biệt về pulse pressure trung bình giữa các nhóm hiệu suất của nam giới
data_male <- subset(data, gender == "M")
# Chạy ANOVA
model_male <- aov(pulse_pressure ~ class, data = data_male)
summary(model_male)
##               Df Sum Sq Mean Sq F value   Pr(>F)    
## class          3   4289  1429.8   12.91 2.05e-08 ***
## Residuals   8453 935831   110.7                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
#P-val = 2.05e-08
# Với mức ý nghĩa alpha = 0.05 , p-val < alpha nên chấp nhận H1 
# Tương tự đối với nhóm nữ : 
#H0 : Không có sự khác biệt về pulse pressure trung bình giữa các nhóm hiệu suất của nữ giới 
#H1 : Có ít nhất 1 sự khác biệt về pulse pressure trung bình giữa các nhóm hiệu suất của nữ giới
data_female <- subset(data, gender == "F")
# Chạy ANOVA
model_female <- aov(pulse_pressure ~ class, data = data_female)
summary(model_female)
##               Df Sum Sq Mean Sq F value Pr(>F)
## class          3    667   222.3   2.036  0.107
## Residuals   4912 536391   109.2

p-val = 0.107 Với mức ý nghĩa alpha =0.05 , p-val >alpha nên chấp nhận H0 Vậy nhận xét ban đầu chưa chính xác và chỉ mang tính ngẫu nhiên vì thực tế , không có sự khác biệt về pulse pressure trung bình giữa các nhóm hiệu suất của nữ giới

6.

ggplot(data, aes(x = bmi_category, y = pulse_pressure, fill = gender)) +
  geom_bar(stat = "summary", fun = "mean", position = "dodge") +
  labs(title = "Pulse pressure (áp lực mạch đập) trung bình theo phân loại BMI và giới tính", 
       x = "Phân loại BMI", y = "Pulse pressure (mmHg)") +
  theme_minimal()

##### Nhận xét : Pulse pressure trung bình ở nam giới của 3 nhóm Normal , Obese , Overweight gần như ngang nhau , của Underweight là thấp nhất, cả 4 nhóm đều cao hơn nữ giới
#Ở nữ giới thì chỉ số lại khá cao trong 2 nhóm Obese và Overweight , nhóm Normal thấp hơn đôi chút và nhóm Underweight là thấp nhất 
# Để kiểm tra nhận xét trên , ta thực hiện các kiểm định sau :
# Ta đã kiểm đinh được Pulse pressure ở nữ giới thấp hơn nam giới trong câu trước , nên giờ ta tiếp tục 

# H0 : Không có sự khác biệt pulse pressure trung bình ở các nhóm BMI ở nam giới 
# H1 : có ít nhất 1 sự khác biệt pulse pressure trung bình ở các nhóm BMI ở nam giới
data_male <- subset(data, gender == "M")
model_male <- aov(pulse_pressure ~ bmi_category, data = data_male)
# Tóm tắt kết quả
summary(model_male)
##                Df Sum Sq Mean Sq F value   Pr(>F)    
## bmi_category    3   2370   790.0   7.121 8.96e-05 ***
## Residuals    8453 937751   110.9                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# P-val = 8.96e-05 
# Với mức ý nghĩa alpha = 0.05 , p-val < alpha nên chấp nhận H1 
# Kế tiếp :
# H0 : Không có sự khác biệt pulse pressure trung bình ở các nhóm BMI ở nữ giới 
# H1 : có ít nhất 1 sự khác biệt pulse pressure trung bình ở các nhóm BMI ở nữ giới

data_female <- subset(data, gender == "F")
model_female <- aov(pulse_pressure ~ bmi_category, data = data_female)
summary(model_female)
##                Df Sum Sq Mean Sq F value Pr(>F)    
## bmi_category    3   8557  2852.2   26.51 <2e-16 ***
## Residuals    4912 528501   107.6                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-val = 2e-16 Với mức ý nghĩa alpha = 0.05 , p-val < alpha nên chấp nhận H1 Vậy nhận xét ban đầu có thể xem là mang ý nghĩa thống kê

3.4.2 .Theo độ tuổi

7. age và pulse_pressure dựa trên class

ggplot(data, aes(x = age, y = pulse_pressure, color = class)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "loess", se = TRUE) +
  labs(title = "Quan hệ giữa Tuổi và Sự khác biệt Pulse pressure theo Class",
       x = "Tuổi", y = "Pulse pressure (mmHg)") +
  scale_color_manual(values = c("A" = "blue", "B" = "red", "C" = "green", "D" = "yellow")) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

##### Nhận xét : dựa vào đường cong thể hiện xu hướng tăng giảm của các chấm dữ liệu , ta thấy chỉ số áp lực mạch đập có xu hướng tăng lên khi tuổi
#tác càng cao 
# Để kiểm định nhận xét trên , ta thực hiện : 
#H0: Không có mối quan hệ tuyến tính giữa tuổi tác và áp lực mạch đập 
#H1: Có mối quan hệ tuyến tính giữa tuổi tác và áp lực mạch đập 
# Chạy mô hình hồi quy
model <- lm(pulse_pressure ~ age, data = data)
# Xem tóm tắt kết quả
summary(model)
## 
## Call:
## lm(formula = pulse_pressure ~ age, data = data)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -44.778  -7.677  -0.271   7.078  86.802 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 47.735704   0.265654  179.69   <2e-16 ***
## age          0.101418   0.006775   14.97   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 10.67 on 13371 degrees of freedom
## Multiple R-squared:  0.01648,    Adjusted R-squared:  0.01641 
## F-statistic: 224.1 on 1 and 13371 DF,  p-value: < 2.2e-16

p-val < 2e-16 Với mức ý nghĩa alpha = 0.05 , p-val < alpha nên ta chấp nhận H1 : Có mối quan hệ tuyến tính giữa tuổi tác và áp lực mạch đập , nghĩa là tuổi tác tăng thì áp lực mạch đập cũng tăng . ta kết luận nhận xét ban đầu là chính xác và mang ý nghĩa thống kê

8. fitness score và age dựa trên class

ggplot(data, aes(x = age, y = fitness_score, color = class)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "loess", se = TRUE) +
  labs(title = "Quan hệ giữa Tuổi và Điểm thể lực theo Class",
       x = "Tuổi", y = "Điểm thể lực") +
  scale_color_manual(values = c("A" = "blue", "B" = "red", "C" = "green", "D" = "yellow")) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

##### Nhận xét : dựa vào các đường cong thể hiện xu hướng tăng giảm của dữ liệu , ta thấy ở cả 4 class thì chỉ số fitness score cao nhất ở giai đoạn
# khoảng 25 tới 35 tuổi . Từ 40 tuổi trở đi thì chỉ số bắt đầu giảm dần. 
# Để kiểm định nhận xét trên , ta thực hiện : 
# Chia ra các nhóm tuổi 
data$age_group <- cut(data$age, 
                      breaks = c(0, 25, 35, 40, Inf), 
                      labels = c("Dưới 25", "25-35", "35-40", "Trên 40"), 
                      right = FALSE)

#H0 : Không có sự khác biệt về fitness score ở các nhóm tuổi 
#H1 :  Có sự khác biệt về fitness score ở các nhóm tuổi
# Tính trung bình fitness_score theo nhóm tuổi
aggregate(fitness_score ~ age_group, data = data, mean)
##   age_group fitness_score
## 1   Dưới 25      89.05074
## 2     25-35      92.15208
## 3     35-40      88.26432
## 4   Trên 40      73.13964
# Kiểm định 
classes <- unique(data$class) 
for (class in classes) {
  cat("\nClass:", class, "\n")
  data_class <- subset(data, class == class)
  
  # Chạy ANOVA
  model <- aov(fitness_score ~ age_group, data = data_class)
  
  # Tóm tắt kết quả
  print(summary(model))
}
## 
## Class: C 
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## age_group       3  959456  319819    1099 <2e-16 ***
## Residuals   13369 3889498     291                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: A 
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## age_group       3  959456  319819    1099 <2e-16 ***
## Residuals   13369 3889498     291                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: B 
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## age_group       3  959456  319819    1099 <2e-16 ***
## Residuals   13369 3889498     291                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: D 
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## age_group       3  959456  319819    1099 <2e-16 ***
## Residuals   13369 3889498     291                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-val <2e-16 Mức ý nghĩa alpha = 0.05 Với các p-val lần lượt : Nhóm A : 2e-16 < alpha Nhóm B : 2e-16 < alpha Nhóm C : 2e-16 < alpha Nhóm D : 2e-16 < alpha -> Ta đủ cơ sở chấp nhận H1 Vậy kết luận nhận xét ban đầu có thể tạm xem là mang ý nghĩa thống kê

9.Quan hệ grip_force và age dựa trên Class

ggplot(data, aes(x = age, y = grip_force, color = class)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "loess", se = TRUE) +
  labs(title = "Quan hệ giữa Tuổi và Lực kẹp theo Class",
       x = "Tuổi", y = "Lực kẹp (kg)") +
  scale_color_manual(values = c("A" = "blue", "B" = "red", "C" = "green", "D" = "yellow")) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

# Nhận xét : Lực kẹp ở cả 4 class đạt đỉnh ở giai đoạn 25 tới 40 tuổi , sau 40 thì chỉ số này có xu hướng giảm xuống 
# Để kiểm tra nhận xét trên có mang tính thống kê hay không , ta có : 
#H0 : Không có sự khác biệt về lực kẹp ở các nhóm tuổi 
#H1 :  Có sự khác biệt về lực kẹp ở các nhóm tuổi
data$age_group <- cut(data$age, 
                      breaks = c(0, 25, 40, Inf), 
                      labels = c("Dưới 25", "25-40", "Trên 40"), 
                      right = FALSE)

# Tính trung bình lực kẹp (gripForce) theo nhóm tuổi và class
aggregate(grip_force ~ age_group + class, data = data, mean)
##    age_group class grip_force
## 1    Dưới 25     A   34.95523
## 2      25-40     A   41.62105
## 3    Trên 40     A   36.40716
## 4    Dưới 25     B   36.12432
## 5      25-40     B   40.84968
## 6    Trên 40     B   35.57088
## 7    Dưới 25     C   37.03683
## 8      25-40     C   39.09221
## 9    Trên 40     C   33.78551
## 10   Dưới 25     D   35.59085
## 11     25-40     D   37.97014
## 12   Trên 40     D   31.83816
classes <- unique(data$class)

for (class in classes) {
  cat("\nClass:", class, "\n")
  data_class <- subset(data, class == class)
  
  # Chạy ANOVA
  model <- aov(grip_force ~ age_group, data = data_class)
  
  # Tóm tắt kết quả
  print(summary(model))
}
## 
## Class: C 
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## age_group       2   92738   46369   439.2 <2e-16 ***
## Residuals   13370 1411435     106                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: A 
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## age_group       2   92738   46369   439.2 <2e-16 ***
## Residuals   13370 1411435     106                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: B 
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## age_group       2   92738   46369   439.2 <2e-16 ***
## Residuals   13370 1411435     106                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: D 
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## age_group       2   92738   46369   439.2 <2e-16 ***
## Residuals   13370 1411435     106                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mức ý nghĩa alpha = 0.05 Với các p-val lần lượt : Nhóm A : 2e-16 < alpha Nhóm B : 2e-16 < alpha Nhóm C : 2e-16 < alpha Nhóm D : 2e-16 < alpha -> Ta đủ cơ sở chấp nhận H1 Vậy nhận xét ban đầu có thể tạm xem là mang ý nghĩa thống kê

10.Quan hệ giữa age và grip_force dựa trên class và gender

ggplot(data, aes(x = age, y = grip_force, color = gender)) +
  geom_point(alpha = 0.7, size = 3) + # Điểm dữ liệu
  geom_smooth(method = "loess", se = TRUE) + # Đường xu hướng
  facet_wrap(~ class) + # Tách biểu đồ theo giới tính
  labs(
    title = "Quan hệ giữa Tuổi và Lực kẹp theo Class và Gender",
    x = "Tuổi",
    y = "Lực kẹp (kg)"
  ) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Nhận xét : Chỉ số lực kẹp ở nam giới nhìn chung tốt hơn nữ giới ở cả 4 class . Tuy nhiên khi tuổi càng cao
thì lực kẹp ở nam giới có xu hướng giảm nhưng ở nữ giới lại không thể hiện xu hướng tăng giảm rõ ràng

Để kiểm tra nhận xét trên có mang ý nghĩa thống kê không , ta thực hiện :
- Dựa vào kết quả kiểm định của plot 9 , ta đã có thể khẳng định chỉ số lực kẹp của nam tốt hơn nữ ở cả 4 class và đồng thời đạt đỉnh ở độ tuổi 25-40 , sau đó thì đi xuống .Tuy nhiên để kiểm tra xu hướng lực kẹp theo tuổi tác ở riêng nữ giới thì ta có : H0 : Không có sự khác biệt về lực kẹp nữ giới ở các nhóm tuổi
H1 : Có ít nhất 1 sự khác biệt về lực kẹp nữ giới ở các nhóm tuổi
data_female <- subset(data, gender == “F”)

# Phân nhóm tuổi (nếu chưa có cột age_group)
data_female$age_group <- cut(data_female$age, 
                             breaks = c(0, 25, 40, Inf), 
                             labels = c("Dưới 25", "25-40", "Trên 40"), 
                             right = FALSE)

# Tính trung bình lực kẹp (grip_force) theo nhóm tuổi và class
aggregate(grip_force ~ age_group + class, data = data_female, mean)
##    age_group class grip_force
## 1    Dưới 25     A   28.33242
## 2      25-40     A   28.99610
## 3    Trên 40     A   26.87683
## 4    Dưới 25     B   26.79236
## 5      25-40     B   27.22451
## 6    Trên 40     B   25.13240
## 7    Dưới 25     C   25.34755
## 8      25-40     C   25.28545
## 9    Trên 40     C   24.13455
## 10   Dưới 25     D   22.93551
## 11     25-40     D   24.10746
## 12   Trên 40     D   23.16845
# Lặp qua từng class để chạy ANOVA
classes <- unique(data_female$class)

for (class in classes) {
  cat("\nClass:", class, "\n")
  data_class <- subset(data_female, class == class)
  
  model <- aov(grip_force ~ age_group, data = data_class)
  print(summary(model))
}
## 
## Class: B 
##               Df Sum Sq Mean Sq F value Pr(>F)    
## age_group      2   5392  2695.9     130 <2e-16 ***
## Residuals   4913 101881    20.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: D 
##               Df Sum Sq Mean Sq F value Pr(>F)    
## age_group      2   5392  2695.9     130 <2e-16 ***
## Residuals   4913 101881    20.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: C 
##               Df Sum Sq Mean Sq F value Pr(>F)    
## age_group      2   5392  2695.9     130 <2e-16 ***
## Residuals   4913 101881    20.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Class: A 
##               Df Sum Sq Mean Sq F value Pr(>F)    
## age_group      2   5392  2695.9     130 <2e-16 ***
## Residuals   4913 101881    20.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Mức ý nghĩa alpha = 0.05
Với các p-val lần lượt :
Nhóm A : 2e-16 < alpha
Nhóm B : 2e-16 < alpha
Nhóm C : 2e-16 < alpha
Nhóm D : 2e-16 < alpha
-> Ta đủ cơ sở chấp nhận H1
> Vậy nhận xét ban đầu có thể tạm xem là mang ý nghĩa thống kê

#11.Quan hệ giữa tuổi và BMI dựa trên giới tính
ggplot(data, aes(x = age, y = bmi, color = gender)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "loess", se = TRUE) +
  labs(
    title = "Quan hệ giữa Tuổi và BMI theo Giới tính",
    x = "Tuổi",
    y = "BMI"
  ) + scale_color_manual(values = c("F" = "blue", "M" = "red")) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

##### Nhận xét : Chỉ số BMI của nam giới ở mức cao trong độ tuổi từ 25 tới 40 . Sau 40 thì chỉ số bắt đầu giảm xuống
# Tuy nhiên ở nữ giới thì chỉ số bmi lại thấp ở độ tuổi 20 tới trước 40 , nhưng ngoài 40 thì lại có xu hướng tăng lên 
# , điều này khá lạ , nguyên nhân có thể do cột Female bị thiếu dữ liệu . 
# Biểu đồ và nhận xét trên chỉ mang tính tham khảo chứ chưa mang ý nghĩa thống kê 

#12.Quan hệ giữa Tuổi và Phần trăm Mỡ Cơ thể theo class
ggplot(data, aes(x = age, y = body_fat, color = gender)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "loess", se = TRUE) +
  facet_wrap(~ class) +
  labs(title = "Quan hệ giữa Tuổi và Body Fat theo Class",
       x = "Tuổi", y = "Body Fat (%)") +
  scale_color_manual(values = c("F" = "blue", "M" = "red")) +
  theme_classic()
## `geom_smooth()` using formula = 'y ~ x'

Nhận xét : % Mỡ(body_fat) ở cả 2 giới đều có xu hướng đi lên khi ở độ tuổi ngoài 40 , ở nữ giới thể hiện xu hướng này rõ hơn

Để kiểm tra nhận xét trên , ta có :
H0 : % mỡ không có sự khác biệt trên cả 2 giới ở mọi độ tuổi
H1 : % mỡ có sự khác biệt từ độ tuổi 40 trở đi trên cả 2 giới

data$age_group <- cut(data$age, 
                      breaks = c(0, 40, Inf), 
                      labels = c("Dưới 40", "Trên 40"), 
                      right = FALSE)

# Tính trung bình % mỡ cơ thể theo giới tính và nhóm tuổi
aggregate(body_fat ~ age_group + gender, data = data, mean)
##   age_group gender body_fat
## 1   Dưới 40      F 27.22851
## 2   Trên 40      F 30.19521
## 3   Dưới 40      M 19.28429
## 4   Trên 40      M 21.95728
data_male <- subset(data, gender == "M")
data_female <- subset(data, gender == "F")

# Kiểm định t cho nam giới
t_test_male <- t.test(body_fat ~ age_group, data = data_male)
cat("Kết quả kiểm định t cho nam giới:\n")
## Kết quả kiểm định t cho nam giới:
print(t_test_male)
## 
##  Welch Two Sample t-test
## 
## data:  body_fat by age_group
## t = -20.565, df = 6226.9, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Dưới 40 and group Trên 40 is not equal to 0
## 95 percent confidence interval:
##  -2.927800 -2.418187
## sample estimates:
## mean in group Dưới 40 mean in group Trên 40 
##              19.28429              21.95728
# Kiểm định t cho nữ giới
t_test_female <- t.test(body_fat ~ age_group, data = data_female)
cat("Kết quả kiểm định t cho nữ giới:\n")
## Kết quả kiểm định t cho nữ giới:
print(t_test_female)
## 
##  Welch Two Sample t-test
## 
## data:  body_fat by age_group
## t = -16.696, df = 4146.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means between group Dưới 40 and group Trên 40 is not equal to 0
## 95 percent confidence interval:
##  -3.315064 -2.618333
## sample estimates:
## mean in group Dưới 40 mean in group Trên 40 
##              27.22851              30.19521

cho alpha =0.05
Cả 2 p-val đều = 2.2e-16 < alpha nên ta chấp nhận H1
> Vậy nhận xét ban đầu có mang ý nghĩa thống kê

14. Quan hệ giữa Systolic và Diastolic theo giới tính

ggplot(data, aes(x = systolic, y = diastolic, color = gender)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Quan hệ giữa Huyết áp Tâm thu và Tâm trương",
       x = "Huyết áp Tâm thu (mmHg)", y = "Huyết áp Tâm trương (mmHg)") +
  theme_light()
## `geom_smooth()` using formula = 'y ~ x'

Nhận xét : Ở cả 2 giới thì huyết áp tâm thu và huyết áp tâm trương đều có xu hướng tỉ lệ thuận , đồng thời 2 chỉ số này ở nam cao hơn nữ

Tách thành 2 biểu đồ nhỏ

14.1 Biểu đồ violin thể hiện phân phối huyết áp tâm thu giữa nam và nữ.

ggplot(data, aes(x = gender, y = systolic, fill = gender)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Gender", y = "Systolic(mmHg)") +
  ggtitle("Phân phối huyết áp tâm thu giữa nam và nữ.") +
  theme_bw()

Qua biểu đồ ta thấy rõ được huyết áp tâm thu của nữ giới thấp hơn nam giới.
Ta sẽ kiểm định “Huyết áp tâm thu của nữ giới thấp hơn nam giới hay không?”
H0: Có sự giống nhau về huyết áp tâm thu giữa nam và nữ.
H1: Huyết áp tâm thu của nữ thấp hơn nam.

diff <- perm_fun(x = data$systolic, nA = 4916, nB = 8457, R = 1000)
mean_a <- mean(data$systolic[data$gender == "F"])
mean_b <- mean(data$systolic[data$gender == "M"])
mean(diff < (mean_a - mean_b))
## [1] 0

p-value = 0 < 0.05 -> Chấp nhận H1: Huyết áp tâm thu của nữ thấp hơn nam.

14.2 Biểu đồ violin thể hiện phân phối huyết áp tâm trương giữa nam và nữ.

ggplot(data, aes(x = gender, y = diastolic, fill = gender)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Gender", y = "Diastolic(mmHg)") +
  ggtitle("Phân phối huyết áp tâm trương giữa nam và nữ.") +
  theme_bw()

Qua biểu đồ ta thấy có vẻ huyết áp tâm trương của nữ giới thấp hơn nam giới.

Ta sẽ kiểm định “Huyết áp tâm trương của nữ giới thấp hơn nam giới hay không?

H0: Có sự giống nhau về huyết áp tâm trương giữa nam và nữ.
H1: Huyết áp tâm trương của nữ thấp hơn nam.

diff <- perm_fun(x = data$diastolic, nA = 4916, nB = 8457, R = 1000)
mean_a <- mean(data$diastolic[data$gender == "F"])
mean_b <- mean(data$diastolic[data$gender == "M"])
mean(diff < (mean_a - mean_b))
## [1] 0
# p-value = 0 < 0.05 -> Chấp nhận H1: Huyết áp tâm trương của nữ thấp hơn nam.

#H0 :  Không có mối quan hệ tuyến tính giữa huyết áp tâm thu và huyết áp tâm trương
#H1 :  có mối quan hệ tuyến tính giữa huyết áp tâm thu và huyết áp tâm trương
cor.test(data$systolic, data$diastolic)
## 
##  Pearson's product-moment correlation
## 
## data:  data$systolic and data$diastolic
## t = 106.87, df = 13371, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6694848 0.6877691
## sample estimates:
##       cor 
## 0.6787321
# Phân tích cho nam
male_data <- filter(data, gender == "M")
cor.test(male_data$systolic, male_data$diastolic)
## 
##  Pearson's product-moment correlation
## 
## data:  male_data$systolic and male_data$diastolic
## t = 75.732, df = 8455, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6228728 0.6482756
## sample estimates:
##       cor 
## 0.6357463
# Phân tích cho nữ
female_data <- filter(data, gender == "F")
cor.test(female_data$systolic, female_data$diastolic)
## 
##  Pearson's product-moment correlation
## 
## data:  female_data$systolic and female_data$diastolic
## t = 66.723, df = 4914, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
##  0.6744902 0.7038352
## sample estimates:
##       cor 
## 0.6894455
# Cả 2 p-val = 2.2e-16 < alpha nên ta chấp nhận H1 
# Vậy kết luận ban đầu mang ý nghĩa thống kê 

15. Quan hệ giữa Cân nặng và áp lực mạch đập theo bmi_category

ggplot(data, aes(x = weight_kg, y = pulse_pressure, color = bmi_category)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "lm", se = TRUE) +
  labs(title = "Quan hệ giữa Cân nặng và Áp lực mạch đập(Pulse pressure) theo BMI",
       x = "Cân nặng (kg)", y = "Áp lực mạch đập (mmHg)") 
## `geom_smooth()` using formula = 'y ~ x'

Nhận xét : Ở cả 4 phân nhóm BMI đều không thể hiện xu hướng rõ ràng cho mối quan hệ giữa cân nặng và áp lực mạch đập . Ta không dám chắc nếu cân năng tặng thì áp lực mạch đập sẽ tăng hoặc ngược lại

15’. Biểu đồ violin thể hiện phân phối áp lực mạch đập theo bmi_category

ggplot(data, aes(x = bmi_category, y = pulse_pressure, fill = bmi_category)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "bmi_category", y = "pulse_pressure (mmHg)") +
  ggtitle("Phân phối áp lực mạch đập giữa các nhóm bmi.") +
  theme_bw()

Qua biểu đồ ta thấy được sự khác biệt về áp lực mạch đập giữa các nhóm bmi
> Ta sẽ kiểm định “Áp lực mạch đập giữa các nhóm bmi có sự khác biệt hay không?”

H0: Không có sự khác biệt về áp lực mạch đập giữa các nhóm bmi.
H1: Có sự khác biệt về áp lực mạch đập giữa các nhóm bmi.

h <- aovp(formula = pulse_pressure ~ bmi_category, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h)
## Component 1 :
##                 Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## bmi_category     3    27081    9027.1 5000 < 2.2e-16 ***
## Residuals    13369  1521623     113.8                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# p-value < 2.2e-16 < 0.05 -> Chấp nhận H1: Có sự khác biệt về áp lực mạch đập giữa các nhóm bmi.

# Biểu đồ tương quan 
corr_matrix <- cor(data[, c("age", "height_cm", "weight_kg","body_fat","diastolic","systolic", 
                            "grip_force","sit_and_bend_forward_cm","sit_ups_counts", "broad_jump_cm", "bmi",
                            "fitness_score","pulse_pressure")], use = "complete.obs")
corrplot(corr_matrix, method = "circle", tl.cex = 0.8)

16.So sánh cân nặng và chiều cao theo giới tính

ggplot(data, aes(x = weight_kg, y = height_cm, color = gender)) +
  geom_point(alpha = 0.7, size = 3) +
  geom_smooth(method = "loess", se = TRUE) +
  labs(
    title = "Quan hệ giữa cân nặng và chiều cao theo Giới tính",
    x = "weight(kg)",
    y = "height(cm)"
  ) + scale_color_manual(values = c("F" = "blue", "M" = "red")) +
  theme_minimal()
## `geom_smooth()` using formula = 'y ~ x'

Nhận xét : Nhìn chung chỉ số cân nặng và chiều cao ở nam giới cao hơn nữ giới . Đường cong thể hiện xu hướng ở nữ giới có dạng gần giống parabol , cho thấy chiều cao tăng lên theo cân nặng lúc ban đầu, sau đó dường như đạt đỉnh nhưng bắt đầu giảm ở các mức cân nặng cao hơn ở nam giới đường cong có xu hướng tăng lên từ từ, cho thấy chiều cao tăng theo cân nặng và có vẻ tiếp tục tăng lên ở mức cân nặng lớn hơn

Xem xét cân nặng giữa nam và nữ

16’. Biểu đồ violin thể hiện phân phối cân nặng giữa nam và nữ

ggplot(data, aes(x = gender, y = weight_kg, fill = gender)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Gender", y = "Cân nặng (kg)") +
  ggtitle("Phân phối cân nặng giữa nam và nữ.") +
  theme_bw()

Qua biểu đồ ta thấy được cân nặng của nữ giới nhẹ hơn nam giới. Ta sẽ kiểm dịnh “Nữ giới có cân nặng nhẹ hơn nam giới

H0: cân nặng giữa nam và nữ bằng nhau.
H1: Cân nặng của nữ thấp hơn nam.

diff <- perm_fun(x = data$weight_kg, nA = 4916, nB = 8457, R = 1000)
mean_a <- mean(data$weight_kg[data$gender == "F"])
mean_b <- mean(data$weight_kg[data$gender == "M"])
mean(diff < (mean_a - mean_b))
## [1] 0

p-value = 0 < 0.05 -> Chấp nhận H1: Cân nặng của nữ thấp hơn nam.
Tuy nhiên chưa tìm được cách kiểm định đường cong parabol có mang ý nghĩa thống kê không ?
Nên nhận xét ban đầu chỉ mang tính tham khảo chứ chưa có ý nghĩa thống kê

17.Box so sánh sit_and_bend_forward_cm theo các nhóm bmi

ggplot(data, aes(x = bmi_category, y = sit_and_bend_forward_cm, fill = bmi_category)) +
  geom_boxplot() +
  theme_minimal() +
  labs(
    title = "So sánh sit_and_bend_forward_cm theo nhóm BMI",
    x = "Nhóm BMI",
    y = "sit_and_bend_forward_cm (cm)"
  )

Nhận xét : Chỉ số sit_and_bend_forward_cm trung bình thấp nhất ở 2 nhóm Obese và Overweight . Nhóm Normal cao hơn và nhóm Underweight là cao nhất

Ta sẽ kiểm định “lực kẹp giữa các nhóm bmi có giống nhau không?”
H0: sit_and_bend_forward_cm giữa các nhóm bmi không khác nhau.
H1: sit_and_bend_forward_cm giữa các nhóm bmi khác nhau.

h <- aovp(formula = grip_force ~ bmi_category, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h)
## Component 1 :
##                 Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## bmi_category     3   173366     57789 5000 < 2.2e-16 ***
## Residuals    13369  1330808       100                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
# p-value < 2.2e-16 < 0.05 -> Chấp nhận H1: sit_and_bend_forward_cm giữa các nhóm bmi khác nhau.

########### theo các class (biểu đồ tổng hợp )
# Normalization (min-max scaling) cho tất cả các cột
data_scaled <- data %>%
  mutate(across(c(bmi, body_fat, fitness_score, pulse_pressure), 
                ~ scales::rescale(.), 
                .names = "{col}_scaled"))

# Vẽ biểu đồ tổng hợp
data_long_scaled <- data_scaled %>%
  gather(key = "variable", value = "value", bmi_scaled, body_fat_scaled, fitness_score_scaled, pulse_pressure_scaled)

ggplot(data_long_scaled, aes(x = class, y = value, fill = variable)) +
  geom_bar(stat = "summary", fun = "mean", position = "dodge") +
  labs(title = "Biểu đồ tổng hợp các chỉ số sau khi scale theo Class", 
       x = "Nhóm Phân loại", 
       y = "Giá trị trung bình (đã chuẩn hóa)") +
  theme_minimal() +
  scale_fill_manual(values = c("bmi_scaled" = "blue", "body_fat_scaled" = "red", 
                               
                               "fitness_score_scaled" = "darkgreen", "pulse_pressure_scaled" = "orange"))

Nhận xét : -Chỉ số BMI ở class A thấp nhất , tuy nhiên ở class B lại cao hơn class C , class D cao nhất
-Fitness score giảm dần đều theo từ class A sang D
-Body fat tăng dần đều từ class A sang D
-Pulse pressure (áp lực mạch đập) gần như ngang nhau ở cả 4 class

18. “Chỉ số nhảy xa của nữ yếu hơn nam hay không?

H0: Chỉ số nhảy xa giữa nam và nữ bằng nhau.
H1: Chỉ số nhảy xa của nữ yếu hơn nam.

# Biểu đồ violin
ggplot(data, aes(x = gender, y = broad_jump_cm, fill = gender)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Gender", y = "Broad Jump(cm)") +
  ggtitle("Phân phối chỉ số nhảy xa giữa nam và nữ.") +
  theme_bw()

Qua biểu đồ, ta thấy sự khác biệt về chỉ số nhảy xa giữa nam và nữ

diff_h6 <- perm_fun(x = data$broad_jump_cm, nA = 4916, nB = 8457, R = 1000)
h6_mean_a <- mean(data$broad_jump_cm[data$gender == "F"])
h6_mean_b <- mean(data$broad_jump_cm[data$gender == "M"])
mean(diff_h6 < (h6_mean_a - h6_mean_b))
## [1] 0

p-value = 0 < 0.05 -> Chấp nhận H1: Chỉ số nhảy xa của nữ yếu hơn nam.

19. Nhóm huyết áp phụ thuộc vào nhóm bmi hay không?

H0: Nhóm huyết áp không phụ thuộc vào nhóm bmi. (2 nhóm độc lập nhau)
H1: Nhóm huyết áp phụ thuộc vào nhóm bmi.

pivot_blood_pressure_bmi <- data %>% 
  count(blood_pressure_category, bmi_category) %>% 
  spread(key = bmi_category, value = n, fill = 0)

# Chuyển đổi bảng pivot thành ma trận
pivot_matrix <- as.matrix(pivot_blood_pressure_bmi[, -1])
rownames(pivot_matrix) <- pivot_blood_pressure_bmi$blood_pressure_category
pivot_matrix
##                      Normal Obese Overweight Underweight
## Elevated               1546    24        456          63
## Hypertension Stage 1   3874   175       1775         104
## Hypertension Stage 2   1415   127        903          19
## Hypotension               5     0          1           0
## Normal                 2391    16        328         151
chisq.test(pivot_matrix)
## Warning in stats::chisq.test(x, y, ...): Chi-squared approximation may be
## incorrect
## 
##  Pearson's Chi-squared test
## 
## data:  pivot_matrix
## X-squared = 803.87, df = 12, p-value < 2.2e-16

p-value < 2.2e-16 < 0.05 -> Chấp nhận H1: Nhóm huyết áp phụ thuộc vào nhóm BMI.

20. “Có hay không sự khác biệt về Chỉ số nhảy bật xa giữa các nhóm bmi.

H0: Không có sự khác biệt về chỉ số nhảy bật xa giữa các nhóm bmi.
H1: Có sự khác biệt về chỉ số nhảy bật xa giữa các nhóm bmi.

# Vẽ biểu đồ violin
ggplot(data, aes(x = bmi_category, y = broad_jump_cm, fill = bmi_category)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "bmi_category", y = "Broad Jump") +
  ggtitle("Phân phối chỉ số nhảy xa giữa các nhóm bmi.") +
  theme_bw()

Qua biểu đồ ta thấy có sự khác biệt về chỉ số nhảy bật xa giữa các nhóm bmi.

h_13 <- aovp(formula = broad_jump_cm ~ bmi_category, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h_13)
## Component 1 :
##                 Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## bmi_category     3   507766    169255 5000 < 2.2e-16 ***
## Residuals    13369 20399757      1526                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value < 2.2e-16 < 0.05 -> Chấp nhận H1: Có sự khác biệt về Chỉ số nhảy bật xa giữa các nhóm bmi

21. “Ở các nhóm chỉ số huyết áp thì độ tuổi có sự khác biệt hay không?

H0: Ở các nhóm chỉ số huyết áp thì độ tuổi không khác biệt
H1: Ở các nhóm chỉ số huyết áp thì độ tuổi có sự khác biệt

# Vẽ biểu đồ violin
ggplot(data, aes(x = blood_pressure_category, y = age, fill = blood_pressure_category)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Blood Pressure Category", y = "Age") +
  ggtitle("Phân phối độ tuổi giữa các nhóm chỉ số huyết áp.") +
  theme_bw()

Qua biểu đồ ta thấy rõ được sự khác biệt về độ tuổi ở các nhóm chỉ số huyết áp

h_14 <- aovp(formula = age ~ blood_pressure_category, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h_14)
## Component 1 :
##                            Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## blood_pressure_category     4    70718   17679.6 5000 < 2.2e-16 ***
## Residuals               13368  2411228     180.4                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value < 2.2e-16 < 0.05 -> Chấp nhận H1: Độ tuổi ở các nhóm chỉ số huyết áp có khác biệt.

22. “Ở các nhóm chỉ số huyết áp thì cân nặng có sự khác biệt hay không?

H0: Cân nặng ở các nhóm chỉ số huyết áp không khác biệt.
H1: Cân nặng ở các nhóm chỉ số huyết áp có khác biệt.

# Vẽ biểu đồ violin
ggplot(data, aes(x = blood_pressure_category, y = weight_kg, fill = blood_pressure_category)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Blood Pressure Category", y = "Weight") +
  ggtitle("Phân phối cân nặng giữa các nhóm chỉ số huyết áp.") +
  theme_bw()

Qua biểu đồ ta thấy rõ được sự khác biệt về cân nặng ở các nhóm chỉ số huyết áp

h_15 <- aovp(formula = weight_kg ~ blood_pressure_category, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h_15)
## Component 1 :
##                            Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## blood_pressure_category     4   201855     50464 5000 < 2.2e-16 ***
## Residuals               13368  1707001       128                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value < 2.2e-16 < 0.05 -> Chấp nhận H1: Cân nặng ở các nhóm chỉ số huyết áp có khác biệt.

23. Ở các nhóm chỉ số huyết áp thì sit_and_bend_forward_cm có sự khác biệt hay không?

H0: sit_and_bend_forward_cm ở các nhóm chỉ số huyết áp không khác biệt.
H1: sit_and_bend_forward_cm ở các nhóm chỉ số huyết áp có khác biệt.

# Vẽ biểu đồ violin
ggplot(data, aes(x = blood_pressure_category, y = sit_and_bend_forward_cm, fill = blood_pressure_category)) +
  geom_violin() +
  geom_boxplot(width = 0.15) +
  labs(x = "Blood Pressure Category", y = "sit_and_bend_forward_cm") +
  ggtitle("Phân phối sit_and_bend_forward_cm giữa các nhóm chỉ số huyết áp.") +
  theme_bw()

Qua biểu đồ ta thấy rõ được sự khác biệt về cân nặng ở các nhóm chỉ số huyết áp

h_15 <- aovp(formula = sit_and_bend_forward_cm ~ blood_pressure_category, data = data, perm = "Prob")
## [1] "Settings:  unique SS "
summary(h_15)
## Component 1 :
##                            Df R Sum Sq R Mean Sq Iter  Pr(Prob)    
## blood_pressure_category     4     7641   1910.18 5000 < 2.2e-16 ***
## Residuals               13368   881038     65.91                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value < 2.2e-16 < 0.05 -> Chấp nhận H1: Sit bend forward ở các nhóm chỉ số huyết áp có khác biệt.

24. Phân lớp hiệu suất phụ thuộc vào giới tính hay không?

H0: Phân lớp hiệu suất không phụ thuộc vào giới tính (2 nhóm độc lập nhau)
H1: Phân lớp hiệu suất phụ thuộc vào giới tính

pivot_class_gender <- data %>%
  count(class, gender) %>%
  spread(key = gender, value = n, fill = 0)

# Chuyển đổi bảng pivot thành ma trận
pivot_matrix <- as.matrix(pivot_class_gender[, -1])
rownames(pivot_matrix) <- pivot_class_gender$class
pivot_matrix
##      F    M
## A 1484 1862
## B 1185 2159
## C 1109 2234
## D 1138 2202
chisq.test(pivot_matrix)
## 
##  Pearson's Chi-squared test
## 
## data:  pivot_matrix
## X-squared = 114.34, df = 3, p-value < 2.2e-16

p-value < 2.2e-16 < 0.05 -> Chấp nhận H1: Phân lớp hiệu suất phụ thuộc vào giới tính

Tổng kết phần 1 : Các biểu đồ và kết quả của những kiểm định phía trên thể hiện những góc nhìn ở nhiều khía cạnh khác nhau trong dữ liệu Điều này có thể góp ích ít nhiều trong việc là những lời khuyên hữu ích cho các chuyên gia sức khỏe .

4 PHẦN 2 : BOOTSTRAP

# Function to calculate the mean for bootstrap
boot_mu_fun <- function(data, ind) {
  data_new <- data[ind]
  out <- mean(data_new)
  return(out)
}

# Sử dụng dữ liệu đã có
data_boot <- data  

# Các biến cần phân tích
variables <- c("age", "height_cm", "weight_kg","body_fat","diastolic","systolic", 
               "grip_force","sit_and_bend_forward_cm","sit_ups_counts", "broad_jump_cm", "bmi",
               "fitness_score","pulse_pressure")

# Khởi tạo bảng kết quả
results <- data.frame(
  Variable = character(0),
  Mean = numeric(0),
  Bias = numeric(0),
  StdError = numeric(0),
  CI_Lower = numeric(0),
  CI_Upper = numeric(0)
)

# Khởi tạo danh sách để lưu đồ thị
plot_list <- list()

for (var in variables) {
  # Thực hiện bootstrap
  out <- boot(data = data_boot[[var]], statistic = boot_mu_fun, R = 1000)
  
  # Trích xuất kết quả
  mean_val <- out$t0
  bias_val <- mean(out$t) - out$t0
  std_error <- sd(out$t)
  
  # Tính khoảng tin cậy
  ci <- boot.ci(out, type = "perc", conf = 0.95)
  ci_lower <- ci$percent[4]  # Giá trị giới hạn dưới
  ci_upper <- ci$percent[5]  # Giá trị giới hạn trên
  
  # Thêm kết quả vào bảng
  results <- rbind(results, data.frame(
    Variable = var,
    Mean = mean_val,
    Bias = bias_val,
    StdError = std_error,
    CI_Lower = ci_lower,
    CI_Upper = ci_upper
  ))
  
  # Tạo biểu đồ histogram
  plot_list[[var]] <- ggplot(data = data.frame(t = out$t), mapping = aes(x = t)) +
    geom_histogram(fill = "gray80", color = "black", bins = 20) +
    geom_vline(xintercept = out$t0, color = "blue", linetype = "dashed") +
    xlab(paste(var)) +
    theme_bw() +
    theme(axis.title.y = element_blank())  # Bỏ nhãn của trục Oy
  
}

# In bảng kết quả
print(results)
##                   Variable      Mean          Bias   StdError  CI_Lower
## 1                      age  36.76939  0.0003727660 0.11831068  36.53761
## 2                height_cm 168.56482 -0.0040849697 0.07353319 168.41482
## 3                weight_kg  67.44905 -0.0044557310 0.10176601  67.24275
## 4                 body_fat  23.23548  0.0037570026 0.06097687  23.11842
## 5                diastolic  78.79868  0.0035596874 0.09194573  78.61982
## 6                 systolic 130.26344 -0.0022300905 0.12688469 130.02090
## 7               grip_force  36.97887 -0.0030030487 0.09167431  36.79339
## 8  sit_and_bend_forward_cm  15.18211 -0.0007937112 0.07300413  15.04404
## 9           sit_ups_counts  39.78307 -0.0049668885 0.12025844  39.55299
## 10           broad_jump_cm 190.26180  0.0164187542 0.34167163 189.63916
## 11                     bmi  23.60526  0.0007975788 0.02569452  23.55690
## 12           fitness_score  84.08543 -0.0013468545 0.16155041  83.76581
## 13          pulse_pressure  51.46476  0.0012238241 0.09565763  51.28055
##     CI_Upper
## 1   37.00520
## 2  168.69340
## 3   67.64535
## 4   23.35563
## 5   78.98886
## 6  130.51347
## 7   37.16003
## 8   15.34046
## 9   40.02783
## 10 190.94048
## 11  23.65731
## 12  84.39899
## 13  51.66508
# Kết hợp các biểu đồ thành lưới
grid.arrange(
  grobs = plot_list,
  ncol = 4,
  top = "Bootstrap Results for Selected Variables"
)

glimpse(data_boot)
## Rows: 13,373
## Columns: 18
## $ age                     <dbl> 27, 25, 31, 32, 28, 36, 42, 33, 54, 28, 42, 57…
## $ gender                  <fct> M, M, M, M, M, F, F, M, M, M, M, F, F, M, M, F…
## $ height_cm               <dbl> 172.3, 165.0, 179.6, 174.5, 173.8, 165.4, 164.…
## $ weight_kg               <dbl> 75.24, 55.80, 78.00, 71.10, 67.70, 55.40, 63.7…
## $ body_fat                <dbl> 21.3, 15.7, 20.1, 18.4, 17.1, 22.0, 32.2, 36.9…
## $ diastolic               <dbl> 80, 77, 92, 76, 70, 64, 72, 84, 85, 81, 63, 69…
## $ systolic                <dbl> 130, 126, 152, 147, 127, 119, 135, 137, 165, 1…
## $ grip_force              <dbl> 54.9, 36.4, 44.8, 41.4, 43.5, 23.8, 22.7, 45.9…
## $ sit_and_bend_forward_cm <dbl> 18.4, 16.3, 12.0, 15.2, 27.1, 21.0, 0.8, 12.3,…
## $ sit_ups_counts          <dbl> 60, 53, 49, 53, 45, 27, 18, 42, 34, 55, 68, 0,…
## $ broad_jump_cm           <dbl> 217, 229, 181, 219, 217, 153, 146, 234, 148, 2…
## $ class                   <fct> C, A, C, B, B, B, D, B, C, B, A, D, C, C, C, A…
## $ bmi                     <dbl> 25.34418, 20.49587, 24.18143, 23.34956, 22.412…
## $ bmi_category            <fct> Overweight, Normal, Normal, Normal, Normal, No…
## $ blood_pressure_category <fct> Hypertension Stage 1, Elevated, Hypertension S…
## $ fitness_score           <dbl> 105.57, 100.82, 87.34, 99.32, 96.15, 63.84, 57…
## $ pulse_pressure          <dbl> 50, 49, 60, 71, 57, 55, 63, 53, 80, 75, 47, 37…
## $ age_group               <fct> Dưới 40, Dưới 40, Dưới 40, Dưới 40, Dưới 40, D…

5 PHẦN 3 : CLASSIFICATION (BÀI TOÁN PHÂN LOẠI)

5.1 .Xử lí cơ bản:

# Lọc bỏ các cột không cần thiết
data <- data[, !names(data) %in% c("bmi_category")]
# Bỏ cột bmi_catgory vì cột trên được tạo từ cột bmi, nếu để có thể gây ra hiện tưởng đa cộng tuyến, làm giảm sự chính xác cho mô hình

# Kiểm tra lại dữ liệu sau khi lọc
glimpse(data)
## Rows: 13,373
## Columns: 17
## $ age                     <dbl> 27, 25, 31, 32, 28, 36, 42, 33, 54, 28, 42, 57…
## $ gender                  <fct> M, M, M, M, M, F, F, M, M, M, M, F, F, M, M, F…
## $ height_cm               <dbl> 172.3, 165.0, 179.6, 174.5, 173.8, 165.4, 164.…
## $ weight_kg               <dbl> 75.24, 55.80, 78.00, 71.10, 67.70, 55.40, 63.7…
## $ body_fat                <dbl> 21.3, 15.7, 20.1, 18.4, 17.1, 22.0, 32.2, 36.9…
## $ diastolic               <dbl> 80, 77, 92, 76, 70, 64, 72, 84, 85, 81, 63, 69…
## $ systolic                <dbl> 130, 126, 152, 147, 127, 119, 135, 137, 165, 1…
## $ grip_force              <dbl> 54.9, 36.4, 44.8, 41.4, 43.5, 23.8, 22.7, 45.9…
## $ sit_and_bend_forward_cm <dbl> 18.4, 16.3, 12.0, 15.2, 27.1, 21.0, 0.8, 12.3,…
## $ sit_ups_counts          <dbl> 60, 53, 49, 53, 45, 27, 18, 42, 34, 55, 68, 0,…
## $ broad_jump_cm           <dbl> 217, 229, 181, 219, 217, 153, 146, 234, 148, 2…
## $ class                   <fct> C, A, C, B, B, B, D, B, C, B, A, D, C, C, C, A…
## $ bmi                     <dbl> 25.34418, 20.49587, 24.18143, 23.34956, 22.412…
## $ blood_pressure_category <fct> Hypertension Stage 1, Elevated, Hypertension S…
## $ fitness_score           <dbl> 105.57, 100.82, 87.34, 99.32, 96.15, 63.84, 57…
## $ pulse_pressure          <dbl> 50, 49, 60, 71, 57, 55, 63, 53, 80, 75, 47, 37…
## $ age_group               <fct> Dưới 40, Dưới 40, Dưới 40, Dưới 40, Dưới 40, D…

5.2 .hàm để đánh giá và nhận xét mô hình

eval_multi_class <- function(x) {
  cc <- sum(diag(x))  # Tổng các giá trị trên đường chéo chính (True Positives)
  sc <- sum(x)        # Tổng tất cả các giá trị trong ma trận nhầm lẫn
  pp <- colSums(x)    # Tổng giá trị theo cột (số mẫu trong mỗi lớp dự đoán)
  tt <- rowSums(x)    # Tổng giá trị theo hàng (số mẫu trong mỗi lớp thực tế)
  
  # Precision, Recall cho mỗi lớp
  prec <- diag(x) / colSums(x)    # Precision = TP / (TP + FP)
  recall <- diag(x) / rowSums(x)  # Recall = TP / (TP + FN)
  
  # Tính Macro Precision, Macro Recall, và Macro F1
  macro_prec <- mean(prec, na.rm = TRUE)  # Trung bình Precision của tất cả các lớp
  macro_recall <- mean(recall, na.rm = TRUE)  # Trung bình Recall của tất cả các lớp
  macro_f1 <- 2 * macro_prec * macro_recall / (macro_prec + macro_recall)  # Tính Macro F1
  
  # Accuracy và Kappa
  acc <- cc / sc  # Accuracy = (TP + TN) / (Tổng số mẫu)
  
  # Tính chỉ số Kappa
  kap <- (cc * sc - sum(pp * tt)) / (sc^2 - sum(pp * tt))
  
  # Trả về kết quả
  return(list(Precision = prec, Recall = recall, Accuracy = acc, Kappa = kap, Macro_F1 = macro_f1))
}

5.3 .K-fold

# Cài đặt số lần gấp (folds) = 5 cho cross-validation
train_control <- trainControl(
  method = "cv", 
  number = 5,  # 5 folds
  summaryFunction = defaultSummary,  
  returnResamp = "all",      
  savePredictions = "all"
)

5.4 .NaiveBayes

5.4.1 .Tính ma trận tương quan vì trước khi sử dụng mô hình NaiveByes phải kiểm tra sự tương quan giữa các biến

numeric_data <- data[sapply(data, is.numeric)]
correlation_matrix <- cor(numeric_data)

# Vẽ biểu đồ nhiệt độ để trực quan hóa ma trận tương quan
corrplot(correlation_matrix, method = "color")

# nhận xét có một số biến có sự tương quan cao -> nên mô hình naiveBayes có thể hoạt động không hiệu quả.
set.seed(123)  # Đảm bảo tính tái lập
nb_model <- NaiveBayes(class ~ ., data = data)
nb_pred <- predict(nb_model, newdata = data)

################ in ra kết quả
# Kết hợp lớp và xác suất của tất cả các lớp vào một bảng
result <- data.frame(nb_pred$posterior)
# Thêm cột cho lớp có xác suất cao nhất
result$Max_Probability <- apply(result, 1, max)
result$Predicted_Class <- apply(result, 1, function(x) names(result)[which.max(x)])
# Thêm cột giá trị thực tế
result$True_Class <- data$class

# Hiển thị kết quả
head(result)
##            A         B          C            D Max_Probability Predicted_Class
## 1 0.59406076 0.3199424 0.08368405 0.0023127460       0.5940608               A
## 2 0.64999161 0.2542972 0.09498089 0.0007303076       0.6499916               A
## 3 0.04939704 0.4601599 0.43691992 0.0535230933       0.4601599               B
## 4 0.31790837 0.4850027 0.19358066 0.0035082619       0.4850027               B
## 5 0.82296375 0.1308662 0.04541421 0.0007558809       0.8229637               A
## 6 0.25423716 0.2558669 0.41563591 0.0742600688       0.4156359               C
##   True_Class
## 1          C
## 2          A
## 3          C
## 4          B
## 5          B
## 6          B
# Đánh giá mô hình
conf_matrix<-table( data$class,nb_pred$class)
eval_multi_class(conf_matrix)
## $Precision
##         A         B         C         D 
## 0.5834747 0.4075417 0.4959748 0.6452768 
## 
## $Recall
##         A         B         C         D 
## 0.7196653 0.3361244 0.4238708 0.7014970 
## 
## $Accuracy
## [1] 0.5452778
## 
## $Kappa
## [1] 0.3936997
## 
## $Macro_F1
## [1] 0.5391089

Nhận xét mô hình - Mô hình có độ chính xác là 54.5%. Đây là một mức độ chính xác tương đối thấp, vì trong correlation matrix ta thấy được một số biến có tương quan cao với nhau nên mô hình yếu
- Kappa : 0.393 cho thấy rằng mức độ tương đồng của mô hình và thực tế là khá yếu
- Mô hình có Precision cao nhất ở lớp D tiếp theo là lớp A, trong khi Precision của lớp B và C còn khá thấp, điều này có thể chỉ ra rằng mô hình gặp khó khăn trong việc phân loại chính xác các đối tượng thuộc lớp B và C.
- Mô hình có Recall cao nhất ở lớp A và D, nhưng lớp B và C bị bỏ sót khá nhiều đối tượng, điều này có thể cho thấy mô hình chưa đủ nhạy để nhận diện đầy đủ các đối tượng của các lớp này
- Macro_F1: 0.54, điều này có nghĩa là mô hình có hiệu suất trung bình khi đánh giá tất cả các lớp, và có sự cân bằng giữa Precision và Recall, tuy nhiên vẫn khá thấp.

5.4.2 .Áp dụng K-fold để kiểm tra thêm lần nữa

nb_model <- train(class ~ ., data = data, method = "naive_bayes", # sử dụng "naive_bayes" trong caret
                  trControl = train_control)
print(nb_model)
## Naive Bayes 
## 
## 13373 samples
##    16 predictor
##     4 classes: 'A', 'B', 'C', 'D' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 10699, 10698, 10698, 10698, 10699 
## Resampling results across tuning parameters:
## 
##   usekernel  Accuracy   Kappa    
##   FALSE      0.3951247  0.1935805
##    TRUE      0.4889698  0.3185828
## 
## Tuning parameter 'laplace' was held constant at a value of 0
## Tuning
##  parameter 'adjust' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were laplace = 0, usekernel = TRUE
##  and adjust = 1.
print(nb_model$resample)  # In ra thông tin chi tiết của từng fold
##     Accuracy     Kappa usekernel laplace adjust Resample
## 1  0.4764398 0.3018055      TRUE       0      1    Fold1
## 2  0.3979058 0.1973984     FALSE       0      1    Fold1
## 3  0.4859813 0.3146123      TRUE       0      1    Fold2
## 4  0.3872897 0.1830697     FALSE       0      1    Fold2
## 5  0.4912150 0.3215104      TRUE       0      1    Fold3
## 6  0.3988785 0.1987761     FALSE       0      1    Fold3
## 7  0.4964486 0.3285688      TRUE       0      1    Fold4
## 8  0.3966355 0.1955266     FALSE       0      1    Fold4
## 9  0.4947644 0.3264172      TRUE       0      1    Fold5
## 10 0.3949140 0.1931317     FALSE       0      1    Fold5

Nhận xét tổng quan về các fold
Accuracy khá thấp , chỉ trong khoảng từ 38% đến 50% -> mô hình dự đoán kém kappa rất thấp, chỉ tầm 19% nếu không dùng kernel đến 32% nếu có dùng kernel và nên sử dụng kernel thì mô hình có vẻ tốt hơn nhưng 31% vẫn là khá thấp nên mức độ tương đồng thực tế thấp

5.4.3 Loại bỏ các cột có tương quan cao hơn ngưỡng

threshold <- 0.7  # Đặt ngưỡng tương quan
high_corr_columns <- findCorrelation(correlation_matrix, cutoff = threshold, names = TRUE)
filtered_data <- numeric_data[, !(names(numeric_data) %in% high_corr_columns)]

### Kiểm tra dữ liệu sau khi lọc
print("Các cột giữ lại sau khi loại bỏ tương quan cao:")
## [1] "Các cột giữ lại sau khi loại bỏ tương quan cao:"
print(names(filtered_data))
## [1] "age"                     "body_fat"               
## [3] "diastolic"               "systolic"               
## [5] "sit_and_bend_forward_cm" "sit_ups_counts"         
## [7] "bmi"                     "pulse_pressure"

5.4.4 Thêm cột ‘class’ vào lại để áp dụng mô hình

filtered_data$class <- data$class
# Áp dụng mô hình Naive Bayes
set.seed(123)
nb_model <- NaiveBayes(class ~ ., data = filtered_data)
nb_pred <- predict(nb_model, newdata = filtered_data)

5.4.5 Đánh giá kết quả

conf_matrix <- table(filtered_data$class, nb_pred$class)
eval_multi_class(conf_matrix)
## $Precision
##         A         B         C         D 
## 0.5839136 0.4150943 0.5167087 0.7283537 
## 
## $Recall
##         A         B         C         D 
## 0.7268380 0.3421053 0.4902782 0.7152695 
## 
## $Accuracy
## [1] 0.5686084
## 
## $Kappa
## [1] 0.4248023
## 
## $Macro_F1
## [1] 0.5647945

mô hình được cải thiện rõ rệt Accuracy 57%, kappa 42%,macro F1 56% nhưng vẫn chưa phải là tốt

5.4.6 dùng k-fold để đánh giá

nb_model <- train(class ~ ., data = filtered_data, method = "naive_bayes", # sử dụng "naive_bayes" trong caret
                  trControl = train_control)
print(nb_model)
## Naive Bayes 
## 
## 13373 samples
##     8 predictor
##     4 classes: 'A', 'B', 'C', 'D' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 10699, 10698, 10698, 10698, 10699 
## Resampling results across tuning parameters:
## 
##   usekernel  Accuracy   Kappa    
##   FALSE      0.5659171  0.4212157
##    TRUE      0.5763852  0.4351680
## 
## Tuning parameter 'laplace' was held constant at a value of 0
## Tuning
##  parameter 'adjust' was held constant at a value of 1
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were laplace = 0, usekernel = TRUE
##  and adjust = 1.
print(nb_model$resample)  # In ra thông tin chi tiết của từng fold
##     Accuracy     Kappa usekernel laplace adjust Resample
## 1  0.5744203 0.4325465      TRUE       0      1    Fold1
## 2  0.5654450 0.4205865     FALSE       0      1    Fold1
## 3  0.5708411 0.4277859      TRUE       0      1    Fold2
## 4  0.5532710 0.4043627     FALSE       0      1    Fold2
## 5  0.5839252 0.4452053      TRUE       0      1    Fold3
## 6  0.5641121 0.4188008     FALSE       0      1    Fold3
## 7  0.5745794 0.4327678      TRUE       0      1    Fold4
## 8  0.5704673 0.4272880     FALSE       0      1    Fold4
## 9  0.5781601 0.4375344      TRUE       0      1    Fold5
## 10 0.5762902 0.4350406     FALSE       0      1    Fold5

nhận xét: mô hình đã được cãi tiến, các hệ số đều cao hơn 57% Accuracy nếu có dùng kernel và 55% nếu không dùng kernel 40% nếu không dùng kernel và 43% nếu dùng kernel cho chỉ số Kappa
–> mô hình đã được cải tiến nhưng chưa tốt

5.5 .LDA

lda_model <- lda(class ~ ., data = data)
lda_model
## Call:
## lda(class ~ ., data = data)
## 
## Prior probabilities of groups:
##         A         B         C         D 
## 0.2502056 0.2500561 0.2499813 0.2497570 
## 
## Group means:
##        age   genderM height_cm weight_kg body_fat diastolic systolic grip_force
## A 35.26689 0.5564854  167.8724  64.41769 20.53987  77.90120 129.3207   38.61122
## B 37.07476 0.6456340  168.5776  66.61314 22.04543  78.66746 130.6675   37.91978
## C 36.70176 0.6682620  169.1646  66.75339 22.63887  78.52468 129.9336   36.59846
## D 38.03653 0.6592814  168.6453  72.01905 27.72457  80.10338 131.1335   34.78231
##   sit_and_bend_forward_cm sit_ups_counts broad_jump_cm      bmi
## A               21.344133       47.84668      202.7325 22.72826
## B               17.411606       42.63170      195.3051 23.31761
## C               14.381923       38.71283      188.8400 23.19409
## D                7.577749       29.92413      174.1425 25.18338
##   blood_pressure_categoryHypertension Stage 1
## A                                   0.4240885
## B                                   0.4458732
## C                                   0.4507927
## D                                   0.4523952
##   blood_pressure_categoryHypertension Stage 2
## A                                   0.1682606
## B                                   0.1851077
## C                                   0.1728986
## D                                   0.2107784
##   blood_pressure_categoryHypotension blood_pressure_categoryNormal
## A                       0.0011954573                     0.2408846
## B                       0.0005980861                     0.2027512
## C                       0.0000000000                     0.2168711
## D                       0.0000000000                     0.2026946
##   fitness_score pulse_pressure age_groupTrên 40
## A      91.54179       51.41949        0.3048416
## B      87.02013       52.00000        0.3687201
## C      83.11666       51.40891        0.3727191
## D      74.64710       51.03015        0.4284431
## 
## Coefficients of linear discriminants:
##                                                       LD1          LD2
## age                                         -0.0570063708  0.024888833
## genderM                                      1.3213447893 -2.770327787
## height_cm                                   -0.0180547185 -0.032903746
## weight_kg                                    0.0589088394 -0.001933907
## body_fat                                     0.0224934286  0.019167847
## diastolic                                    0.0030445946  0.010086589
## systolic                                     0.0003717598  0.003845359
## grip_force                                  -0.0381590900  0.025160049
## sit_and_bend_forward_cm                     -0.0980982631 -0.038445146
## sit_ups_counts                              -0.0641404933  0.016933718
## broad_jump_cm                               -0.0010028560  0.011143689
## bmi                                         -0.0403969145  0.274291228
## blood_pressure_categoryHypertension Stage 1  0.0114371565 -0.171538800
## blood_pressure_categoryHypertension Stage 2 -0.0110896169 -0.062795505
## blood_pressure_categoryHypotension          -0.2702223818  3.563473381
## blood_pressure_categoryNormal                0.0159905038  0.330102810
## fitness_score                               -0.0179500992  0.020911715
## pulse_pressure                              -0.0023075844 -0.002829457
## age_groupTrên 40                             0.3315820219 -0.265427714
##                                                       LD3
## age                                          4.878326e-02
## genderM                                      2.192846e+00
## height_cm                                    1.069156e-01
## weight_kg                                   -1.712391e-01
## body_fat                                     5.831815e-02
## diastolic                                   -1.135742e-02
## systolic                                     6.815659e-05
## grip_force                                  -1.486474e-02
## sit_and_bend_forward_cm                      5.105295e-02
## sit_ups_counts                               9.884280e-03
## broad_jump_cm                               -1.351608e-03
## bmi                                          5.408281e-01
## blood_pressure_categoryHypertension Stage 1 -2.066363e-01
## blood_pressure_categoryHypertension Stage 2 -1.305600e-01
## blood_pressure_categoryHypotension           2.158208e-01
## blood_pressure_categoryNormal               -6.488169e-01
## fitness_score                               -1.378863e-03
## pulse_pressure                               1.128867e-02
## age_groupTrên 40                            -5.664235e-01
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.9784 0.0196 0.0020
# Dự đoán trên toàn bộ dữ liệu
lda_pred <- predict(lda_model, newdata = data)
# 
conf_matrix_lda <- table(data$class, lda_pred$class)
# Đánh giá mô hình
eval_multi_class(conf_matrix_lda)
## $Precision
##         A         B         C         D 
## 0.6867639 0.4595420 0.5287709 0.8238482 
## 
## $Recall
##         A         B         C         D 
## 0.7319187 0.4500598 0.5662579 0.7281437 
## 
## $Accuracy
## [1] 0.6190832
## 
## $Kappa
## [1] 0.4921043
## 
## $Macro_F1
## [1] 0.6219004

nhận thấy lda tốt hơn naiveBayes vì nó giả sử các feature có cùng phương sai và dùng phân phối chuẩn để tính ước lượng hàm mật độ xác suất đồng thời fj(x) –> ít bị ảnh hưởng bởi vấn đề độc lập tuyến tính hơn so với Naive Bayes

  • Accuracy 62% cho thấy mô hình dự đoán độ chính xác trung bình
  • Kappa: 49% cho thấy sự tương đồng trung bình
  • F1 62% : tốt hơn Naive Bayes nhưng vẫn chưa thực sự hiệu quả

5.5.1 dùng k-fold để đánh giá

set.seed(123)
lda_model <- train(class ~ ., data = data, 
                   method = "lda",  # Sử dụng method "lda"  trong caret
                   trControl = train_control)

# In kết quả
print(lda_model)
## Linear Discriminant Analysis 
## 
## 13373 samples
##    16 predictor
##     4 classes: 'A', 'B', 'C', 'D' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 10699, 10698, 10698, 10698, 10699 
## Resampling results:
## 
##   Accuracy   Kappa    
##   0.6143726  0.4858251
print(lda_model$resample)  # In ra thông tin chi tiết của từng fold
##    Accuracy     Kappa parameter Resample
## 1 0.6062079 0.4749400      none    Fold1
## 2 0.6123364 0.4831081      none    Fold2
## 3 0.6097196 0.4796280      none    Fold3
## 4 0.6157009 0.4875941      none    Fold4
## 5 0.6278983 0.5038555      none    Fold5

5.6 .multinominal logistic

5.6.1 . Xây dựng mô hình Multinomial Logistic Regression

logistic_model <- multinom(class ~ ., data = data)
## # weights:  84 (60 variable)
## initial  value 18538.914491 
## iter  10 value 14520.430572
## iter  20 value 14303.052390
## iter  30 value 14254.209653
## iter  40 value 13141.143634
## iter  50 value 11732.668643
## iter  60 value 11476.231377
## iter  70 value 11465.902058
## final  value 11465.847203 
## converged
summary(logistic_model)
## Call:
## multinom(formula = class ~ ., data = data)
## 
## Coefficients:
##   (Intercept)         age  genderM   height_cm weight_kg    body_fat
## B    23.36239 -0.09207194 2.578242 -0.07080774 0.1556105 0.004432550
## C    45.92846 -0.17191438 4.277800 -0.15013998 0.3167667 0.005867987
## D    49.43424 -0.25455227 5.579639 -0.16935870 0.3809205 0.067846186
##     diastolic      systolic  grip_force sit_and_bend_forward_cm sit_ups_counts
## B 0.001287344 -0.0009768004 -0.07147075              -0.1674633    -0.07993911
## C 0.003484409 -0.0013018916 -0.11708437              -0.2801331    -0.14716012
## D 0.012453792  0.0033302678 -0.15784162              -0.4357350    -0.23799187
##   broad_jump_cm        bmi blood_pressure_categoryHypertension Stage 1
## B -0.0013300386 -0.2305837                                  0.03797596
## C -0.0006912729 -0.6088756                                  0.09588131
## D  0.0099302550 -0.4892788                                  0.01589422
##   blood_pressure_categoryHypertension Stage 2
## B                                  0.04542059
## C                                  0.05690739
## D                                 -0.06367781
##   blood_pressure_categoryHypotension blood_pressure_categoryNormal
## B                         0.01504006                   -0.10995861
## C                        -8.07588715                   -0.07020669
## D                        -7.55334999                    0.08822504
##   fitness_score pulse_pressure age_groupTrên 40
## B   -0.05381588   -0.002264140        0.4603493
## C   -0.09419674   -0.004786308        0.9902046
## D   -0.13957016   -0.009123523        1.4077356
## 
## Std. Errors:
##    (Intercept)         age   genderM   height_cm   weight_kg    body_fat
## B 0.0009493272 0.005676575 0.1519114 0.004111372 0.008762230 0.007722575
## C 0.0012829111 0.006650624 0.1774031 0.004746087 0.009962578 0.009033720
## D 0.0015520106 0.008236594 0.2179112 0.005691502 0.012054444 0.011474123
##     diastolic    systolic  grip_force sit_and_bend_forward_cm sit_ups_counts
## B 0.002731703 0.002552323 0.005978847             0.006774328    0.003965063
## C 0.003186084 0.002967897 0.006923178             0.007797089    0.004701255
## D 0.004095057 0.003732386 0.008530598             0.009191355    0.005861282
##   broad_jump_cm        bmi blood_pressure_categoryHypertension Stage 1
## B   0.001886693 0.02794967                                  0.09633386
## C   0.002128815 0.03216175                                  0.11220600
## D   0.002572108 0.03914184                                  0.14489247
##   blood_pressure_categoryHypertension Stage 2
## B                                   0.1372163
## C                                   0.1595585
## D                                   0.2039875
##   blood_pressure_categoryHypotension blood_pressure_categoryNormal
## B                       4.495435e-03                     0.1024530
## C                       1.288655e-06                     0.1199064
## D                       6.393069e-07                     0.1553121
##   fitness_score pulse_pressure age_groupTrên 40
## B   0.001996385    0.001999745        0.1338251
## C   0.002390123    0.002318838        0.1551601
## D   0.002919002    0.002941265        0.1937905
## 
## Residual Deviance: 22931.69 
## AIC: 23039.69
# Dự đoán trên toàn bộ dữ liệu
logistic_pred <- predict(logistic_model, newdata = data, type = "class")

5.6.2 Ma trận nhầm lẫn

conf_matrix_logistic <- table(data$class, logistic_pred)
eval_multi_class(conf_matrix_logistic)
## $Precision
##         A         B         C         D 
## 0.6948344 0.4685839 0.5298575 0.7799577 
## 
## $Recall
##         A         B         C         D 
## 0.7396892 0.4482656 0.5228836 0.7736527 
## 
## $Accuracy
## [1] 0.6211022
## 
## $Kappa
## [1] 0.4948007
## 
## $Macro_F1
## [1] 0.6197124

5.6.3 .K fold

logistic_model <- train(class ~ ., data = data, 
                        method = "multinom",  # Sử dụng method "multinom" cho Multinomial Logistic trong caret
                        trControl = train_control)
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11438.067937
## iter  20 value 11363.865668
## iter  30 value 11330.916786
## iter  40 value 10839.213946
## iter  50 value 9522.432600
## iter  60 value 9204.598074
## iter  70 value 9194.106796
## final  value 9193.829702 
## converged
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11438.081295
## iter  20 value 11363.883036
## iter  30 value 11330.975558
## iter  40 value 10837.213595
## iter  50 value 9524.212496
## iter  60 value 9220.543574
## final  value 9220.533369 
## converged
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11438.067950
## iter  20 value 11363.865686
## iter  30 value 11330.916845
## iter  40 value 10839.211985
## iter  50 value 9522.434701
## iter  60 value 9204.681679
## iter  70 value 9194.587600
## final  value 9194.332974 
## converged
## # weights:  84 (60 variable)
## initial  value 14831.963370 
## iter  10 value 11558.011766
## iter  20 value 11422.848160
## iter  30 value 11380.330199
## iter  40 value 10284.935478
## iter  50 value 9460.964938
## iter  60 value 9210.764947
## iter  70 value 9200.769414
## final  value 9200.635614 
## converged
## # weights:  84 (60 variable)
## initial  value 14831.963370 
## iter  10 value 11558.028224
## iter  20 value 11422.869882
## iter  30 value 11380.402880
## iter  40 value 10280.650542
## iter  50 value 9459.462832
## iter  60 value 9226.900670
## final  value 9226.899737 
## converged
## # weights:  84 (60 variable)
## initial  value 14831.963370 
## iter  10 value 11558.011782
## iter  20 value 11422.848182
## iter  30 value 11380.330271
## iter  40 value 10284.930922
## iter  50 value 9460.962842
## iter  60 value 9210.858268
## iter  70 value 9201.262634
## final  value 9201.139631 
## converged
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11406.930552
## iter  20 value 11336.055785
## iter  30 value 11299.449846
## iter  40 value 10786.264242
## iter  50 value 9433.075719
## iter  60 value 9194.189887
## iter  70 value 9186.471346
## final  value 9186.370964 
## converged
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11406.943834
## iter  20 value 11336.073047
## iter  30 value 11299.504728
## iter  40 value 10784.387743
## iter  50 value 9432.622835
## iter  60 value 9208.345907
## final  value 9208.340307 
## converged
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11406.930565
## iter  20 value 11336.055802
## iter  30 value 11299.449901
## iter  40 value 10786.262419
## iter  50 value 9433.075842
## iter  60 value 9194.256133
## iter  70 value 9186.837020
## final  value 9186.745447 
## converged
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11406.411174
## iter  20 value 11332.693073
## iter  30 value 11294.318002
## iter  40 value 10733.100330
## iter  50 value 9413.386485
## iter  60 value 9122.964957
## iter  70 value 9109.532751
## final  value 9109.387712 
## converged
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11406.424743
## iter  20 value 11332.710473
## iter  30 value 11294.371084
## iter  40 value 10729.861402
## iter  50 value 9419.489690
## iter  60 value 9142.160147
## final  value 9142.149161 
## converged
## # weights:  84 (60 variable)
## initial  value 14830.577075 
## iter  10 value 11406.411187
## iter  20 value 11332.693090
## iter  30 value 11294.318054
## iter  40 value 10733.097116
## iter  50 value 9413.393752
## iter  60 value 9123.081211
## iter  70 value 9110.187637
## final  value 9110.054377 
## converged
## # weights:  84 (60 variable)
## initial  value 14831.963370 
## iter  10 value 11390.252252
## iter  20 value 11322.995581
## iter  30 value 11279.900730
## iter  40 value 10772.930941
## iter  50 value 9439.224759
## iter  60 value 9159.117507
## iter  70 value 9147.677129
## final  value 9147.544748 
## converged
## # weights:  84 (60 variable)
## initial  value 14831.963370 
## iter  10 value 11390.266087
## iter  20 value 11323.012766
## iter  30 value 11279.980187
## iter  40 value 10770.883624
## iter  50 value 9439.730158
## iter  60 value 9175.324262
## final  value 9175.315383 
## converged
## # weights:  84 (60 variable)
## initial  value 14831.963370 
## iter  10 value 11390.252266
## iter  20 value 11322.995598
## iter  30 value 11279.900809
## iter  40 value 10772.928930
## iter  50 value 9439.225420
## iter  60 value 9159.212084
## iter  70 value 9148.214411
## final  value 9148.093238 
## converged
## # weights:  84 (60 variable)
## initial  value 18538.914491 
## iter  10 value 14520.455092
## iter  20 value 14303.083684
## iter  30 value 14254.341044
## iter  40 value 13132.451001
## iter  50 value 11734.870952
## iter  60 value 11497.614377
## final  value 11497.611696 
## converged
# In kết quả
print(logistic_model)
## Penalized Multinomial Regression 
## 
## 13373 samples
##    16 predictor
##     4 classes: 'A', 'B', 'C', 'D' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 10698, 10699, 10698, 10698, 10699 
## Resampling results across tuning parameters:
## 
##   decay  Accuracy   Kappa    
##   0e+00  0.6172890  0.4897159
##   1e-04  0.6172142  0.4896162
##   1e-01  0.6177374  0.4903137
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was decay = 0.1.
print(logistic_model$resample)  # In ra thông tin chi tiết của từng fold
##     Accuracy     Kappa decay Resample
## 1  0.6082243 0.4776302 0e+00    Fold1
## 2  0.6112150 0.4816175 1e-01    Fold1
## 3  0.6082243 0.4776302 1e-04    Fold1
## 4  0.6185490 0.4913966 0e+00    Fold2
## 5  0.6159312 0.4879064 1e-01    Fold2
## 6  0.6185490 0.4913965 1e-04    Fold2
## 7  0.6310280 0.5080366 0e+00    Fold3
## 8  0.6269159 0.5025536 1e-01    Fold3
## 9  0.6306542 0.5075381 1e-04    Fold3
## 10 0.6074766 0.4766310 0e+00    Fold4
## 11 0.6127103 0.4836082 1e-01    Fold4
## 12 0.6071028 0.4761324 1e-04    Fold4
## 13 0.6211668 0.4948853 0e+00    Fold5
## 14 0.6219147 0.4958826 1e-01    Fold5
## 15 0.6215408 0.4953840 1e-04    Fold5

nhận xét : - Accuracy Kappa lần lượt à 0.6140564 0.4854008 giống như trên. Các fold có Accuracy và kappa không quá chênh lệch nhau

Sau khi train và sử dụng K-fold thì thấy mô hình phù hợp nhất cho dữ liệu là lqd và mutinomial logistic.

5.7 .các mô hình không có trong bài học

5.7.1 .Decision Tree

# Xây dựng mô hình Decision Tree
tree_model <- rpart(class ~ ., data = data, method = "class")
rpart.plot(tree_model)

# Dự đoán trên toàn bộ dữ liệu
tree_pred <- predict(tree_model, newdata = data, type = "class")

# Ma trận nhầm lẫn
conf_matrix_tree <- table(data$class, tree_pred)

# Đánh giá mô hình
eval_multi_class(conf_matrix_tree)
## $Precision
##         A         B         C         D 
## 0.5377734 0.4154037 0.6872753 0.7919371 
## 
## $Recall
##         A         B         C         D 
## 0.8084280 0.5000000 0.2859707 0.6940120 
## 
## $Accuracy
## [1] 0.5721229
## 
## $Kappa
## [1] 0.4294675
## 
## $Macro_F1
## [1] 0.5895511

nhận xét mô hình Decision Tree có độ chính xác thấp cho bộ dữ liệu này.

5.7.2 .XGBoost

# Chuyển đổi dữ liệu thành dạng ma trận sparse
sparse_matrix <- sparse.model.matrix(class ~ . - 1, data = data)
label_vector <- as.numeric(data$class) - 1  # XGBoost yêu cầu nhãn dạng số từ 0 trở lên

# Xây dựng tập dữ liệu DMatrix
dtrain <- xgb.DMatrix(data = sparse_matrix, label = label_vector)

# Thiết lập tham số mô hình
params <- list(
  objective = "multi:softmax",
  num_class = length(unique(data$class)),  # Số lượng lớp
  eval_metric = "mlogloss",
  max_depth = 6,
  eta = 0.3
)

# Huấn luyện mô hình XGBoost
xgb_model <- xgb.train(
  params = params,
  data = dtrain,
  nrounds = 100
)

# Dự đoán trên toàn bộ dữ liệu
xgb_pred <- predict(xgb_model, newdata = dtrain)
xgb_pred <- as.factor(xgb_pred)  # Chuyển về dạng factor để so sánh
# Ma trận nhầm lẫn
conf_matrix_xgb <- table(data$class, xgb_pred)

# Đánh giá mô hình
eval_multi_class(conf_matrix_xgb)
## $Precision
##         0         1         2         3 
## 0.9090409 0.9296498 0.9710510 0.9993835 
## 
## $Recall
##         A         B         C         D 
## 0.9886432 0.9207536 0.9231229 0.9706587 
## 
## $Accuracy
## [1] 0.9507964
## 
## $Kappa
## [1] 0.9343946
## 
## $Macro_F1
## [1] 0.9515374

Nhận thấy XGboost có độ chính xác rất cao, lên đến 95% Accuracy, 93% Kappa, 94% Marco_F1 cho thấy các chỉ số đều tốt. về Precision và recall thì ở Class B có thấp hơn so với các class còn lại một chút nhưng vẫn là rất tốt.

5.7.3 .Random Forest

# Xây dựng mô hình Random Forest
set.seed(123)  # Đặt seed để kết quả có thể tái lập
rf_model <- randomForest(class ~ ., data = data, ntree = 100, mtry = sqrt(ncol(data) - 1), importance = TRUE)

# Xem thông tin mô hình
print(rf_model)
## 
## Call:
##  randomForest(formula = class ~ ., data = data, ntree = 100, mtry = sqrt(ncol(data) -      1), importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 100
## No. of variables tried at each split: 4
## 
##         OOB estimate of  error rate: 26.85%
## Confusion matrix:
##      A    B    C    D class.error
## A 2804  477   58    7   0.1619845
## B  746 1954  505  139   0.4156699
## C  258  647 2291  147   0.3146874
## D   36  186  385 2733   0.1817365
# Dự đoán trên toàn bộ dữ liệu
rf_pred <- predict(rf_model, newdata = data)

# Ma trận nhầm lẫn
conf_matrix_rf <- table(data$class, rf_pred)

# Đánh giá mô hình
eval_multi_class(conf_matrix_rf)
## $Precision
## A B C D 
## 1 1 1 1 
## 
## $Recall
## A B C D 
## 1 1 1 1 
## 
## $Accuracy
## [1] 1
## 
## $Kappa
## [1] 1
## 
## $Macro_F1
## [1] 1
# Tầm quan trọng của các biến
importance(rf_model)
##                                  A          B         C          D
## age                      37.315415 30.3839846 27.629774  5.9232753
## gender                   18.917236 10.1912145 11.585000  6.3327246
## height_cm                16.581506  9.3446170  9.234518  6.7651478
## weight_kg                22.412603 11.2381619 14.166254 12.4356513
## body_fat                 20.099449 11.9620501 24.139529 23.4815745
## diastolic                 5.373951  1.3990092  2.413622  1.6691067
## systolic                  6.862318  3.7283090  5.824963  1.9709846
## grip_force               31.323368 18.7685906 11.305121  5.1469791
## sit_and_bend_forward_cm 140.418272 56.7242993 55.881633 57.2310271
## sit_ups_counts           66.218884 34.2835955 29.861043 27.3303606
## broad_jump_cm            22.022419  9.5006800 10.359584  5.7378870
## bmi                      20.062422  9.9473180 15.702351 22.2935071
## blood_pressure_category   3.173675  1.5534362  2.274099  0.8254121
## fitness_score            30.506404 14.5552986 15.808395 14.4226870
## pulse_pressure            4.956436 -0.1488105  3.404134 -1.2463281
## age_group                13.406551 12.7730466 11.456509 -0.5615428
##                         MeanDecreaseAccuracy MeanDecreaseGini
## age                                41.704178        731.37782
## gender                             18.968745        130.94041
## height_cm                          20.155443        478.15879
## weight_kg                          26.370342        603.04766
## body_fat                           31.762987        749.72351
## diastolic                           5.868007        347.52231
## systolic                           10.351415        358.03808
## grip_force                         31.681382        612.63930
## sit_and_bend_forward_cm            97.862801       2342.19404
## sit_ups_counts                     53.923931       1079.35795
## broad_jump_cm                      23.642642        497.03326
## bmi                                28.629383        694.54238
## blood_pressure_category             4.003126        130.45906
## fitness_score                      31.699551        824.71679
## pulse_pressure                      4.367276        358.86375
## age_group                          16.194828         90.30791
varImpPlot(rf_model)

# Random forest chính xác nhất ( dựa trên đánh giá lần train này vì các chỉ số đều 1 hết)
# Kiểm tra lại bằng K-fold

rf_model <- train(class ~ ., data = data, 
                  method = "rf",  
                  trControl = train_control)

# In kết quả
print(rf_model)
## Random Forest 
## 
## 13373 samples
##    16 predictor
##     4 classes: 'A', 'B', 'C', 'D' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold) 
## Summary of sample sizes: 10697, 10699, 10698, 10699, 10699 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    2    0.7082198  0.6109554
##   10    0.7459812  0.6613011
##   19    0.7420935  0.6561170
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 10.
print(rf_model$resample)  # In ra thông tin chi tiết của từng fold
##     Accuracy     Kappa mtry Resample
## 1  0.6988042 0.5983986    2    Fold1
## 2  0.7473842 0.6631689   10    Fold1
## 3  0.7380419 0.6507115   19    Fold1
## 4  0.7064323 0.6085723    2    Fold2
## 5  0.7427076 0.6569385   10    Fold2
## 6  0.7363500 0.6484623   19    Fold2
## 7  0.7031776 0.6042324    2    Fold3
## 8  0.7368224 0.6490908   10    Fold3
## 9  0.7334579 0.6446044   19    Fold3
## 10 0.7157816 0.6210370    2    Fold4
## 11 0.7460733 0.6614224   10    Fold4
## 12 0.7445774 0.6594271   19    Fold4
## 13 0.7169035 0.6225368    2    Fold5
## 14 0.7569185 0.6758848   10    Fold5
## 15 0.7580404 0.6773795   19    Fold5

sau khi thực hiện K-fold thì thấy Accuray chỉ khoảng 74% và kappa ở tầm 66% khá tốt, thấp hơn lần trước có thể là do sự ngẫu nhiên hoặc mô hình trước bị overfit
Các chỉ số của các fold đều chênh lệch nhau không quá nhiều. Sử dụng mtry =18 hoặc =10 có vẻ là sự lựa chọn tốt.

Nhận xét : Mô hình tốt nhất là random forest vì độ chính xác khá ổn và tính giải thích được của nó khá tốt. Ngoài ra nhờ vào nó ta có thể xem xét được tầm quan trọng của các biến đối với độ chính xác mô hình.
- Mean Decrease Accuracy (MDA):
Thước đo này đánh giá mức độ giảm độ chính xác của mô hình nếu loại bỏ một biến cụ thể. Biến nào có giá trị MDA cao, tức là nếu loại bỏ biến đó, mô hình sẽ giảm độ chính xác đáng kể, cho thấy biến này rất quan trọng. - Mean Decrease Gini (MDG):
Thước đo này thể hiện mức giảm độ thuần nhất của các node trong cây quyết định khi chia tách dữ liệu theo biến cụ thể. Giá trị MDG cao cho thấy biến đó rất hữu ích để phân tách dữ liệu, và do đó đóng vai trò lớn trong dự đoán.

Dựa vào 2 biểu đồ MDA và MDG , ta xếp tầm quan trọng của các yếu tố (biến) tới hiệu quả việc tập thể dục theo các nhóm như sau :
* Nhóm yếu tố quan trọng cao nhất:
- sit_and_bend_forward_cm:
Trên cả hai thước đo, yếu tố này đều có giá trị cao nhất, vượt trội so với các biến khác. Điều này có thể do tính linh hoạt (flexibility) của cơ thể là một chỉ số quan trọng trong việc đánh giá sức khỏe và hiệu quả tập luyện.
Trong thực tế, khả năng gập người thường phản ánh rõ mức độ dẻo dai và khả năng vận động của cơ thể.
* Nhóm yếu tố quan trọng cao:
- sit_ups_counts:
Đây là bài tập cơ bản giúp kiểm tra sức mạnh cơ bụng, một phần quan trọng trong hiệu quả tập luyện. Việc biến này có giá trị cao ở cả hai thước đo cho thấy nó là một yếu tố then chốt trong việc phân loại hiệu suất.
- age:
Tuổi tác thường liên quan chặt chẽ đến sức khỏe thể chất. Người trẻ thường có thể lực tốt hơn, trong khi tuổi cao có thể ảnh hưởng đến sức bền, sự linh hoạt và hiệu quả luyện tập.
- fitness_score:
Đây có thể là biến tổng hợp hoặc điểm số phản ánh nhiều khía cạnh của thể chất. Một yếu tố tổng quát như thế này có giá trị cao là hợp lý.
* Nhóm yếu tố quan trọng trung bình:
- body_fat, bmi, grip_force, weight_kg, height_cm, broad_jump_cm:
Các yếu tố này có giá trị MDA và MDG trung bình. Tuy nhiên, chúng vẫn quan trọng vì đều có liên quan trực tiếp đến khả năng vận động, sức mạnh và sức bền của cơ thể.
Ví dụ: grip_force phản ánh sức mạnh cơ tay, còn broad_jump_cm đánh giá sức mạnh và sự bùng nổ của cơ chân.
* Nhóm yếu tố quan trọng thấp:
- systolic, diastolic, pulse_pressure, blood_pressure_category, gender:
Các yếu tố này liên quan đến huyết áp và giới tính, tuy có tác động đến sức khỏe tổng thể nhưng ít trực tiếp hơn đối với hiệu quả tập luyện cụ thể.
Ví dụ, huyết áp (systolic/diastolic) quan trọng hơn khi đánh giá nguy cơ bệnh tim mạch hơn là đánh giá khả năng tập luyện. Tương tự, gender thường chỉ ảnh hưởng gián tiếp thông qua các yếu tố khác như hormone hoặc sự khác biệt về cơ sinh học.